; Reads the best storm tracks of extratropical cyclones ; from track.F program for verificaion dataset and GFS output ; and computes an average track and SLP error ; 95% confidence intervals based on paired t-test are made to examine significance ; ; Andrew Kren ; May 11, 2018 begin ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; begin user's settings ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; CTL_name = "CTL_1" ; name and path name for your control, CASE SENSITIVE EXP_name = (/"CTL_1_noIR"/) ; name and path for your experiment, CASE SENSITIVE ; modify paths accordingly for location of track files path_ver = "/scratch4/BMC/qosap/Andrew.Kren/tracker/best/" ; path of verification tracks path_ctl = "/scratch4/BMC/qosap/Andrew.Kren/tracker/" + CTL_name + "/" ; path of CTL tracks path_exper = "/scratch4/BMC/qosap/Andrew.Kren/tracker/" + EXP_name + "/" ; path of EXP tracks fcst_length = 168 ; how far out in hours you want the error plots to go fcst_interval = 12 ; forecast interval, i.e. every 6hrs, 12hrs, etc. start_date = "2016041500" ; start date of Verification dataset end_date = "2016053000" ; end date of Verification dataset, should be longer than GFS end date for correct verification stats end_date_gfs = "2016051500" ; end date of GFS initialization init_interval = 24 ; interval of your GFS initializations, every 6, 12, 24, etc. output_type = "png" ; "x11", "ps", "eps", "pdf", etc. ; format of the output images you want colors_plot = (/"Black","darkorange1","Green4"/) ; colors for plots system("rm -f *." + output_type) ; remove pngs before running program ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; end of user's settings ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; function closest_val_AnyOrder(xVal[1]:numeric, x:numeric) local xAbsDif, xMinVal, iClose begin xAbsDif = abs(xVal-x) iClose = minind(xAbsDif) return(iClose) ; original doc says "first occurrence" end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; begin of main program ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Verification track file file_ver = path_ver + "track.grads.file" ;;; dates to loop over to read GFS tracks for CTL and EXP year_st = stringtointeger(str_get_cols(start_date,0,3)) year_end = year_st + 1 dates_temp = yyyymmddhh_time(year_st,year_end,fcst_interval,"string") ; string dates indice_start = ind(dates_temp.eq.(start_date)) indice_end = ind(dates_temp.eq.(end_date)) dates_loop = stringtointeger(dates_temp(indice_start:indice_end)) dates_loop@units = "YYYYMMDDHH" delete([/dates_temp,indice_start,indice_end/]) dates_temp = yyyymmddhh_time(year_st,year_end,init_interval,"string") ; string dates indice_start = ind(dates_temp.eq.(start_date)) indice_end = ind(dates_temp.eq.(end_date_gfs)) dates_loop_gfs = stringtointeger(dates_temp(indice_start:indice_end)) dates_loop_gfs@units = "YYYYMMDDHH" delete([/dates_temp,indice_start,indice_end/]) ;;; ; read in the best track file print("Reading best tracks...") data1 = asciiread(file_ver,-1,"string") ; space separated values delim = " " storm_num = tointeger(str_get_field(data1(4::),1,delim)) ; skip first four lines due to header time_indice = tointeger(str_get_field(data1(4::),2,delim)) ; indices that correspond to time lon = tofloat(str_get_field(data1(4::),3,delim)) ; lon from -180 to 180 lat = tofloat(str_get_field(data1(4::),4,delim)) ; latitude slp = tofloat(str_get_field(data1(4::),5,delim)) ; SLP in hPa iflag = tointeger(str_get_field(data1(4::),6,delim)) ; iflag that denotes whether storm ongoing (0), genesis (-1), lysis (1), or occurs for only one time period (-9) ; determine number of extratropical storms num_storms = num(storm_num.eq.1.and.iflag.ne.-9) ; each time it resets to 1 its a new storm that is tracked, filter out storms that occur for one time period only indices_start = ind(iflag.eq.-1.and.iflag.ne.-9) ; genesis stage (-1) denotes start of storm tracking indices_end = ind(iflag.ne.-9.and.iflag.eq.1) ; lysis stage (1) denotes end of storm tracking max_length = max(indices_end - indices_start) ; maximum lifetime of storm tracked in time (hrs) to be used for determining array sizes below max_hr = fcst_interval * max_length ; maxiumum hrs for tracking a storm from all storms in file ; define arrays fcst_hr = ispan(0,max_hr,fcst_interval) ; forecast hours for program slp_array = new((/num_storms,dimsizes(fcst_hr)/),"float") lat_array = slp_array lon_array = slp_array lat_lon_array = new((/num_storms,dimsizes(fcst_hr),dimsizes(fcst_hr)/),"float") ymd_array = new((/num_storms,dimsizes(fcst_hr)/),"string") storm_name = ispan(0,num_storms,1) ; storm names as integers from 1 to num_storms ; loop over storms in verification dataset, storing data do i=0,num_storms-1 start_ind = indices_start(i) end_ind = indices_end(i) end_ind_array = end_ind - start_ind ; for correct input into array temp_dates = time_indice(start_ind:end_ind) - 1 ; for array indexing ymd_array(i,0:end_ind_array) = dates_loop(temp_dates) slp_array(i,0:end_ind_array) = slp(start_ind:end_ind) lat_array(i,0:end_ind_array) = lat(start_ind:end_ind) lon_array(i,0:end_ind_array) = lon(start_ind:end_ind) delete(temp_dates) end do delete([/data1,storm_num,time_indice,lon,lat,slp,iflag,indices_start,indices_end/]) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; read in CTL files and compute errors at each initialization ; loop over GFS dates print("Analyzing CTL Tracks...") do i=0,dimsizes(dates_loop_gfs)-1 print("Analyzing initialization date: " + dates_loop_gfs(i)) file_ctl = path_ctl + dates_loop_gfs(i) + "/" + "track.grads.file" data1 = asciiread(file_ctl,-1,"string") ; space separated values delim = " " storm_num = tointeger(str_get_field(data1(4::),1,delim)) ; skip first four lines due to header time_indice = tointeger(str_get_field(data1(4::),2,delim)) ; indices that correspond to time lon = tofloat(str_get_field(data1(4::),3,delim)) ; lon from -180 to 180 lat = tofloat(str_get_field(data1(4::),4,delim)) ; latitude slp = tofloat(str_get_field(data1(4::),5,delim)) ; SLP in hPa iflag = tointeger(str_get_field(data1(4::),6,delim)) ; iflag that denotes whether storm ongoing (0), genesis (-1), lysis (1), or occurs for only one time period (-9) num_storms_ctl = num(storm_num.eq.1.and.iflag.ne.-9) ; each time it resets to 1 its a new storm that is tracked, filter out storms that occur for one time period only indices_start = ind(iflag.eq.-1.and.iflag.ne.-9) ; genesis stage (-1) denotes start of storm tracking indices_end = ind(iflag.ne.-9.and.iflag.eq.1) ; lysis stage (1) denotes end of storm tracking max_length = max(indices_end - indices_start) ; maximum lifetime of storm tracked in time (hrs) to be used for determining array sizes below max_hr = fcst_interval * max_length ; maxiumum hrs for tracking a storm from all storms in file fcst_hr_cycle = ispan(0,fcst_length,fcst_interval) ; forecast hours for program dates_temp = yyyymmddhh_time(year_st,year_end,fcst_interval,"string") ; string dates indice_start = ind(dates_temp.eq.(dates_loop_gfs(i))) dates_cycle = dates_temp(indice_start:indice_start+dimsizes(fcst_hr_cycle)-1) ; define arrays slp_array_ctl = new((/num_storms,dimsizes(fcst_hr)/),"float") lat_array_ctl = slp_array_ctl lon_array_ctl = slp_array_ctl ymd_array_ctl = new((/num_storms,dimsizes(fcst_hr)/),"string") storm_name_ctl = ispan(0,num_storms,1) ; storm names as integers from 1 to num_storms ; loop over storms in ctl file, storing data do j=0,num_storms_ctl-1 start_ind = indices_start(j) end_ind = indices_end(j) end_ind_array = end_ind - start_ind ; for correct input into array temp_dates = time_indice(start_ind:end_ind) - 1 ; for array indexing ymd_array_ctl(j,0:end_ind_array) = dates_cycle(temp_dates) slp_array_ctl(j,0:end_ind_array) = slp(start_ind:end_ind) lat_array_ctl(j,0:end_ind_array) = lat(start_ind:end_ind) lon_array_ctl(j,0:end_ind_array) = lon(start_ind:end_ind) delete(temp_dates) end do ; end loop over file ; match ctl storms to verification storm dataset and compute errors do q=0,num_storms-1 ;TRY COORDINATE SUBSCRIPTING TO FIND CLOSEST VALUE TO VERIFICATION DATASET lon_diff = min(abs(lon_array_ctl(q,0) - lon_array(:,0))) test2 = closest_val_AnyOrder(lon_array_ctl(q,0), lon_array(:,0)) print(test2) test3 = closest_val_AnyOrder(lat_array_ctl(q,0), lat_array(:,0)) print(test3) exit print(ind(lat_array_ctl(q,0) - lat_array(:,0).eq.lat_diff)) exit lon_diff = min(abs(lon_array_ctl(q,0) - lon_array(:,0))) print(lon_diff) exit slp_diff = abs(slp_array_ctl(q,0) - slp_array(:,0)) match_indice = ind(lat_diff.lt.5.and.lon_diff.lt.5.and.slp_diff.lt.10) print(match_indice) delete([/match_indice,lat_diff,lon_diff,slp_diff/]) end do ; end loop over ctl storms exit delete([/data1,storm_num,time_indice,lon,lat,slp,iflag,indices_start,indices_end,fcst_hr_cycle,dates_temp,indice_start,dates_cycle/]) end do ; end loop over gfs init. dates ;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; end of reading file section ;;;;;;;;;;;;;;;;;;;;;;;;;;;;; count_storms = 0 ; determine number of forecasts per ctl and exp for correct results fcst_dates_ctl = get_unique_values(ctl_ymd) total_fcst_dates_ctl = count_unique_values(ctl_ymd) fcst_dates_exp = new((/dimsizes(EXP_name),100/),"string") total_fcst_dates_exp = new((/dimsizes(EXP_name)/),"float") do a=0,dimsizes(EXP_name)-1 temp = get_unique_values(exp_ymd(a,:,:)) fcst_dates_exp(a,0:dimsizes(temp)-1) = temp total_fcst_dates_exp(a) = count_unique_values(exp_ymd(a,:,:)) delete(temp) end do ; compute storm track and intensity error over each storm individually and save for computing statistics over all storms, including t-test print("Computing track and intensity errors for CTL...") ; first for CTL only since the dates between CTL and EXP are not always the same so treat separately do i=0,dimsizes(g5nr_name)-1 ; first make sure the G5NR name is also found in the CTL experiment, if not, skip the storm storm_index = ind(ctl_name.eq.g5nr_name(i)) if (.not.ismissing(storm_index)) then ; then storm exists in both NR and CTL count_storms = count_storms + 1 print("Analyzing storm: " + g5nr_name(i) + "...") ; find first instance of matching dates between starting NR track and starting GFS run date, if exists fcst_dates = get_unique_values(ctl_ymd(storm_index,:)) ; fcst init dates from GFS fcst_runs = count_unique_values(ctl_ymd(storm_index,:)) ; number of GFS runs for this particular TC storm do j=0,fcst_runs-1 ; loop over forecast runs check = ind(str_squeeze(g5nr_ymd(i,:)).eq.str_squeeze(fcst_dates(j))) if (.not.ismissing(check)) then start_date = str_squeeze(g5nr_ymd(i,check)) ; starting date for NR start_index = check ; index start of nature run print("Analyzing GFS run " + start_date + "...") ; loop over time periods of NR and GFS, computing and saving storm track and intensity error fcst_indices = ind(str_squeeze(ctl_ymd(storm_index,:)).eq.start_date) ; find all indices where the GFS run is initialized at the NR starting date if (.not.all(ismissing(fcst_indices))) then ; loop over fcsts of GFS, computing storm track and intensity error from NR ; conform length of NR tracks or GFS run to each other to account for runs to persist longer or shorter than the best track (if GFS forecast dies/goes beyond best track) length_gfs = dimsizes(fcst_indices) ; number of forecast hours in this GFS run length_nr = num(.not.ismissing(g5nr_ymd(i,start_index::))) ; number of dates in NR from start of GFS run to the end of best track dataset for this storm if (length_gfs.gt.length_nr) then ; means there are more gfs fcst hrs than best tracks start_ind_nr = start_index start_ind_gfs = fcst_indices(0) end_ind_nr = num(.not.ismissing(g5nr_ymd(i,::)))-1 ; last date index of G5NR best track file end_ind_gfs = (end_ind_nr - start_ind_nr) + fcst_indices(0) ; last date index of G5NR best track file end if if (length_nr.ge.length_gfs) then ; means there are more nr tracks than gfs fcst hrs start_ind_nr = start_index start_ind_gfs = fcst_indices(0) end_ind_gfs = fcst_indices(num(.not.ismissing(fcst_indices))-1) ; last date index of GFS track file end_ind_nr = (end_ind_gfs - start_ind_gfs) + start_index ; last date index of G5NR best track file end if ; track and intensity error end_2 = end_ind_gfs - start_ind_gfs ;print(ctl_lat(storm_index,start_ind_gfs:end_ind_gfs)) ;print(ctl_lon(storm_index,start_ind_gfs:end_ind_gfs)) ;print(g5nr_lat(i,start_ind_nr:end_ind_nr)) ;print(g5nr_lon(i,start_ind_nr:end_ind_nr)) ctl_ymd_array(i,j) = ctl_ymd(storm_index,start_ind_gfs) track_error_ctl(i,j,0:end_2) = abs(gc_latlon(ctl_lat(storm_index,start_ind_gfs:end_ind_gfs),ctl_lon(storm_index,start_ind_gfs:end_ind_gfs),g5nr_lat(i,start_ind_nr:end_ind_nr),g5nr_lon(i,start_ind_nr:end_ind_nr),2,4)) ; returns km distance intensity_error_ctl(i,j,0:end_2) = abs( ctl_pres(storm_index,start_ind_gfs:end_ind_gfs) - g5nr_pres(i,start_ind_nr:end_ind_nr) ) wind_error_ctl(i,j,0:end_2) = abs( ctl_wind(storm_index,start_ind_gfs:end_ind_gfs) - g5nr_wind(i,start_ind_nr:end_ind_nr) ) ctl_lat_array(i,j,0:end_2) = ctl_lat(storm_index,start_ind_gfs:end_ind_gfs) ctl_lon_array(i,j,0:end_2) = ctl_lon(storm_index,start_ind_gfs:end_ind_gfs) delete(fcst_indices) end if ; end if for checking missing values of fcst_indices else ; then no matching date with NR and GFS for this init. cycle print("No matching NR date with GFS Init Run of " + fcst_dates(j) + "...") end if ; end if for checking matching dates end do ; end loop over GFS forecast runs delete(fcst_dates) delete(fcst_runs) else ; else for checking if both storms exist in NR and CTL print("Skipping storm: " + g5nr_name(i) + "...") end if ; end if for checking if both storms exist in NR and CTL end do ; end loop over computing stats over each storm of NR ;-------------------------------------------------------------------------------------------------------- print("Computing track and intensity errors for EXP(s)...") ; now for EXP only since the dates between CTL and EXP are not always the same so treat separately do a=0,dimsizes(EXP_name)-1 print("Analyzing EXP: " + EXP_name(a) + "...") do i=0,dimsizes(g5nr_name)-1 ; first make sure the G5NR name is also found in the CTL experiment, if not, skip the storm storm_index = ind(exp_name(a,:).eq.g5nr_name(i)) if (.not.ismissing(storm_index)) then ; then storm exists in both NR and EXP print("Analyzing storm: " + g5nr_name(i) + "...") ; find first instance of matching dates between starting NR track and starting GFS run date, if exists fcst_dates = get_unique_values(exp_ymd(a,storm_index,:)) ; fcst init dates from GFS fcst_runs = count_unique_values(exp_ymd(a,storm_index,:)) ; number of GFS runs for this particular TC storm do j=0,fcst_runs-1 ; loop over forecast runs check = ind(str_squeeze(g5nr_ymd(i,:)).eq.str_squeeze(fcst_dates(j))) if (.not.ismissing(check)) then start_date = str_squeeze(g5nr_ymd(i,check)) ; starting date for NR start_index = check ; index start of nature run print("Analyzing GFS run " + start_date + "...") ; loop over time periods of NR and GFS, computing and saving storm track and intensity error fcst_indices = ind(str_squeeze(exp_ymd(a,storm_index,:)).eq.start_date) ; find all indices where the GFS run is initialized at the NR starting date if (.not.all(ismissing(fcst_indices))) then ; loop over fcsts of GFS, computing storm track and intensity error from NR ; conform length of NR tracks or GFS run to each other to account for runs to persist longer or shorter than the best track (if GFS forecast dies/goes beyond best track) length_gfs = dimsizes(fcst_indices) ; number of forecast hours in this GFS run length_nr = num(.not.ismissing(g5nr_ymd(i,start_index::))) ; number of dates in NR from start of GFS run to the end of best track dataset for this storm if (length_gfs.gt.length_nr) then ; means there are more gfs runs than best tracks start_ind_nr = start_index start_ind_gfs = fcst_indices(0) end_ind_nr = num(.not.ismissing(g5nr_ymd(i,::)))-1 ; last date index of G5NR best track file end_ind_gfs = (end_ind_nr - start_ind_nr) + fcst_indices(0) ; last date index of G5NR best track file end if if (length_nr.ge.length_gfs) then ; means there are more nr tracks than gfs runs start_ind_nr = start_index start_ind_gfs = fcst_indices(0) end_ind_gfs = fcst_indices(num(.not.ismissing(fcst_indices))-1) ; last date index of GFS track file end_ind_nr = (end_ind_gfs - start_ind_gfs) + start_index ; last date index of G5NR best track file end if ; track and intensity error end_2 = end_ind_gfs - start_ind_gfs ;if (g5nr_name(i).eq."AL02") then ;print(exp_lat(storm_index,start_ind_gfs:end_ind_gfs)) ;print(exp_lon(storm_index,start_ind_gfs:end_ind_gfs)) ;print(g5nr_lat(i,start_ind_nr:end_ind_nr)) ;print(g5nr_lon(i,start_ind_nr:end_ind_nr)) ;end if exp_ymd_array(a,i,j) = exp_ymd(a,storm_index,start_ind_gfs) track_error_exp(a,i,j,0:end_2) = abs(gc_latlon(exp_lat(a,storm_index,start_ind_gfs:end_ind_gfs),exp_lon(a,storm_index,start_ind_gfs:end_ind_gfs),g5nr_lat(i,start_ind_nr:end_ind_nr),g5nr_lon(i,start_ind_nr:end_ind_nr),2,4)) ; returns km distance intensity_error_exp(a,i,j,0:end_2) = abs( exp_pres(a,storm_index,start_ind_gfs:end_ind_gfs) - g5nr_pres(i,start_ind_nr:end_ind_nr) ) wind_error_exp(a,i,j,0:end_2) = abs( exp_wind(a,storm_index,start_ind_gfs:end_ind_gfs) - g5nr_wind(i,start_ind_nr:end_ind_nr) ) exp_lat_array(a,i,j,0:end_2) = exp_lat(a,storm_index,start_ind_gfs:end_ind_gfs) exp_lon_array(a,i,j,0:end_2) = exp_lon(a,storm_index,start_ind_gfs:end_ind_gfs) delete(fcst_indices) end if ; end if for checking missing values of fcst_indices else ; then no matching date with NR and GFS for this init. cycle print("No matching NR date with GFS Init Run of " + fcst_dates(j) + "...") end if ; end if for checking matching dates end do ; end loop over GFS forecast runs delete(fcst_dates) delete(fcst_runs) else ; else for checking if both storms exist in NR and EXP print("Skipping storm: " + g5nr_name(i) + "...") end if ; end if for checking if both storms exist in NR and EXP end do ; end loop over computing stats over each storm of NR end do ; end loop over each EXP name ; save data to netcdf if warranted if (create_netcdf.eq."YES") then filo = "saved_track_files" + ".nc" system("rm -rf "+filo) ; remove any pre-existing file qf = addfile(filo,"c") qf@title = "Saved Track Files" qf@source = filo qf@creation_date = systemfunc("date") print("Creating saved netcdf file...") qf->g5nr_lat = g5nr_lat qf->g5nr_lon = g5nr_lon qf->ctl_lat = ctl_lat(0:1,0:28) qf->ctl_lon = ctl_lon(0:1,0:28) qf->exp_lat = exp_lat(:,0:1,0:28) qf->exp_lon = exp_lon(:,0:1,0:28) end if ; ;;;;;;;;;;;;;;;;;;;;;;;; ;; plots results for each storm and GFS forecast separately ;;;;;;;;;;;;;;;;;;;;;;;; print("Plotting results as function of Storm and GFS run...") check_exp = new((/dimsizes(EXP_name)/),"integer") ; now for EXP only since the dates between CTL and EXP are not always the same so treat separately do i=0,dimsizes(g5nr_name)-1 ; first make sure the G5NR name is also found in the CTL experiment, if not, skip the storm storm_index = ind(exp_name(0,:).eq.g5nr_name(i)) if (.not.ismissing(storm_index)) then ; then storm exists in both NR and EXPs print("plotting storm: " + g5nr_name(i) + "...") ; find first instance of matching dates between starting NR track and starting GFS run dates, if exists all_fcsts = array_append_record(fcst_dates_ctl,ndtooned(fcst_dates_exp),0) all_uniq = get_unique_values(all_fcsts) all_count = count_unique_values(all_fcsts) fcst_dates = all_uniq ; fcst init dates from GFS fcst_runs = all_count ; number of forecast init dates from all CTL and EXP runs regardless of storm do j=0,fcst_runs-1 ; loop over forecast runs check = ind(str_squeeze(g5nr_ymd(i,:)).eq.str_squeeze(fcst_dates(j))) check_ctl = ind(str_squeeze(ctl_ymd_array(i,:)).eq.str_squeeze(fcst_dates(j))) do a=0,dimsizes(EXP_name)-1 check_exp(a) = ind(str_squeeze(exp_ymd_array(a,i,:)).eq.str_squeeze(fcst_dates(j))) end do if (.not.ismissing(check) .and. .not.ismissing(check_ctl) .and. .not.any(ismissing(check_exp))) then ; results must exist for all cases, G5NR, CTL, and EXP(s), or nothing will be plotted, account for varying gfs fcsts start_date = str_squeeze(g5nr_ymd(i,check)) ; starting date for NR start_index = check ; index start of nature run if (.not.all(ismissing(track_error_ctl(i,check_ctl,:))) .and. .not.any(ismissing(check_exp)) ) then ; then make track and intensity error plots if (num(.not.ismissing(track_error_ctl(i,check_ctl,:))).gt.1) then ; need 2 elements to plot, just check one exp case for now print("Plotting GFS run " + start_date + "...") ; count number of non-missing values, conform each to the smaller array count_ctl = num(.not.ismissing(track_error_ctl(i,check_ctl,:))) count_exp1 = new((/dimsizes(EXP_name)/),"integer") do a=0,dimsizes(EXP_name)-1 if (.not.ismissing(check_exp(a))) then count_exp1(a) = num(.not.ismissing(track_error_exp(a,i,check_exp(a),:))) end if end do count_g5 = num(.not.ismissing(g5nr_lon(i,start_index::))) count_exp = max(count_exp1) ; find max to plot the full range of all, confined by CTL, of course! if (count_ctl.le.count_exp) then x_axis_max = tointeger(count_ctl) - 1 ; subtract one based on array indexing x_axis_max_temp = tointeger(x_axis_max) + 1 end if if (count_exp.lt.count_ctl) then x_axis_max = tointeger(count_exp) - 1 ; subtract one based on array indexing x_axis_max_temp = tointeger(x_axis_max) + 1 end if ; sizes for y-axes bounds max1 = max(track_error_ctl(i,check_ctl,:)) max2 = max(track_error_exp(:,i,check_exp,:)) maxes = (/max1,max2/) max_track = max(maxes) max1 = max(wind_error_ctl(i,check_ctl,:)) max2 = max(wind_error_exp(:,i,check_exp,:)) maxes = (/max1,max2/) max_wind = max(maxes) max1 = max(intensity_error_ctl(i,check_ctl,:)) max2 = max(intensity_error_exp(:,i,check_exp,:)) maxes = (/max1,max2/) max_intensity = max(maxes) ; when plotting tracks on the map, make sure the length of all are the same for better comparison if (x_axis_max_temp.gt.count_g5) then length_plot_ctl = count_g5 - 1 length_plot_g5 = start_index + (count_g5 - 1) end if if (count_g5.ge.x_axis_max_temp) then length_plot_g5 = start_index + (count_ctl - 1) length_plot_ctl = count_ctl - 1 end if ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; plot errors as a function of forecast lead time ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; wks = gsn_open_wks (output_type,"track_error_fcst_hour-" + g5nr_name(i) + "_" + start_date) res = True ; plot mods desired res@tiMainString = "Track Error (" + g5nr_name(i) + "): " + start_date res@trYMaxF = max_track + 100.0 res@trYMinF = 0 res@gsnFrame = False ; don't advance frame yet ;colors_plot = (/"Black","Red","Orange","Green","Blue","Yellow"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@vpWidthF = 0.8 ; set width and height res@vpHeightF = 0.3 res@vpXF = 0.15 res@vpYF = 0.9 res@trXMinF = 0 res@trXMaxF = x_axis_max x_axis = ispan(0,tointeger(x_axis_max),1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour res@tmXBValues = ispan(0,tointeger(x_axis_max),2) res@tmXBLabels = fcst_hr(0:tointeger(x_axis_max):2) res@tiYAxisString = "Track Error (km)" res@tiXAxisString = "Forecast lead time (hr)" res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.012 res@tmYLLabelFontHeightF = 0.015 res@tmXMajorGrid = True ; implement x grid res@tmXMajorGridThicknessF = 1.0 ; 2.0 is default res@tmXMajorGridLineDashPattern = 2 res@xyLineThicknessF = 8.0 res@tiMainFontHeightF = 0.021 res@xyLineColor = colors_plot(0) plot0 = gsn_csm_xy(wks,x_axis,track_error_ctl(i,check_ctl,0:x_axis_max),res) ; plot additional tracks res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 8.0 ;print(x_axis) ;print(count_exp1(0)-1) if (count_exp1(0).gt.count_ctl) then plot1 = gsn_csm_xy(wks,x_axis,track_error_exp(0,i,check_exp(0),0:x_axis_max),res) else plot1 = gsn_csm_xy(wks,x_axis(0:count_exp1(0)-1),track_error_exp(0,i,check_exp(0),0:count_exp1(0)-1),res) end if overlay(plot0,plot1) if (num_exp_cases.ge.2) then res@xyLineColor = colors_plot(2) res@xyLineThicknessF = 8.0 if (count_exp1(1).gt.count_ctl) then plot2 = gsn_csm_xy(wks,x_axis,track_error_exp(1,i,check_exp(1),0:x_axis_max),res) else plot2 = gsn_csm_xy(wks,x_axis(0:count_exp1(1)-1),track_error_exp(1,i,check_exp(1),0:count_exp1(1)-1),res) end if overlay(plot0,plot2) end if if (num_exp_cases.ge.3) then res@xyLineColor = colors_plot(3) res@xyLineThicknessF = 8.0 if (count_exp1(2).gt.count_ctl) then plot3 = gsn_csm_xy(wks,x_axis,track_error_exp(2,i,check_exp(2),0:x_axis_max),res) else plot3 = gsn_csm_xy(wks,x_axis(0:count_exp1(2)-1),track_error_exp(2,i,check_exp(2),0:count_exp1(2)-1),res) end if overlay(plot0,plot3) end if if (num_exp_cases.ge.4) then res@xyLineColor = colors_plot(4) res@xyLineThicknessF = 8.0 if (count_exp1(3).gt.count_ctl) then plot4 = gsn_csm_xy(wks,x_axis,track_error_exp(3,i,check_exp(3),0:x_axis_max),res) else plot4 = gsn_csm_xy(wks,x_axis(0:count_exp1(3)-1),track_error_exp(3,i,check_exp(3),0:count_exp1(3)-1),res) end if overlay(plot0,plot4) end if xtxt2 = (/0.19,0.19,0.19,0.19,0.19/) ytxt2 = (/0.64,0.62,0.6,0.58,0.56/) labels2 = array_append_record(CTL_name,EXP_name,0) ;(/CTL_name,EXP_name/) ; Attach a legend txres = True txres@txFontHeightF = 0.016 do k=0,dimsizes(labels2)-1 txres@txFontColor = colors_plot(k) gsn_text_ndc(wks,labels2(k),xtxt2(k),ytxt2(k),txres) end do draw(plot0) frame(wks) delete(res) wks = gsn_open_wks (output_type,"wind_error_fcst_hour-" + g5nr_name(i) + "_" + start_date) res = True ; plot mods desired res@tiMainString = "Wind Error (" + g5nr_name(i) + "): " + start_date res@trYMaxF = max_wind + 10.0 res@trYMinF = 0 res@gsnFrame = False ; don't advance frame yet ;colors_plot = (/"Black","Red"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@vpWidthF = 0.8 ; set width and height res@vpHeightF = 0.3 res@vpXF = 0.15 res@vpYF = 0.9 res@trXMinF = 0 res@trXMaxF = x_axis_max x_axis = ispan(0,tointeger(x_axis_max),1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour res@tmXBValues = ispan(0,tointeger(x_axis_max),2) res@tmXBLabels = fcst_hr(0:tointeger(x_axis_max):2) res@tiYAxisString = "Wind Error (kts)" res@tiXAxisString = "Forecast lead time (hr)" res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.012 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 8.0 res@tiMainFontHeightF = 0.021 res@xyLineColor = colors_plot(0) res@tmXMajorGrid = True ; implement x grid res@tmXMajorGridThicknessF = 1.0 ; 2.0 is default res@tmXMajorGridLineDashPattern = 2 plot0 = gsn_csm_xy(wks,x_axis,wind_error_ctl(i,check_ctl,0:x_axis_max),res) ; plot additional tracks res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 8.0 if (count_exp1(0).gt.count_ctl) then plot1 = gsn_csm_xy(wks,x_axis,wind_error_exp(0,i,check_exp(0),0:x_axis_max),res) else plot1 = gsn_csm_xy(wks,x_axis(0:count_exp1(0)-1),wind_error_exp(0,i,check_exp(0),0:count_exp1(0)-1),res) end if overlay(plot0,plot1) if (num_exp_cases.ge.2) then res@xyLineColor = colors_plot(2) res@xyLineThicknessF = 8.0 if (count_exp1(1).gt.count_ctl) then plot2 = gsn_csm_xy(wks,x_axis,wind_error_exp(1,i,check_exp(1),0:x_axis_max),res) else plot2 = gsn_csm_xy(wks,x_axis(0:count_exp1(1)-1),wind_error_exp(1,i,check_exp(1),0:count_exp1(1)-1),res) end if overlay(plot0,plot2) end if if (num_exp_cases.ge.3) then res@xyLineColor = colors_plot(3) res@xyLineThicknessF = 8.0 if (count_exp1(2).gt.count_ctl) then plot3 = gsn_csm_xy(wks,x_axis,wind_error_exp(2,i,check_exp(2),0:x_axis_max),res) else plot3 = gsn_csm_xy(wks,x_axis(0:count_exp1(2)-1),wind_error_exp(2,i,check_exp(2),0:count_exp1(2)-1),res) end if overlay(plot0,plot3) end if if (num_exp_cases.ge.4) then res@xyLineColor = colors_plot(4) res@xyLineThicknessF = 8.0 if (count_exp1(3).gt.count_ctl) then plot4 = gsn_csm_xy(wks,x_axis,wind_error_exp(3,i,check_exp(3),0:x_axis_max),res) else plot4 = gsn_csm_xy(wks,x_axis(0:count_exp1(3)-1),wind_error_exp(3,i,check_exp(3),0:count_exp1(3)-1),res) end if overlay(plot0,plot4) end if ; Attach a legend txres = True txres@txFontHeightF = 0.016 do k=0,dimsizes(labels2)-1 txres@txFontColor = colors_plot(k) gsn_text_ndc(wks,labels2(k),xtxt2(k),ytxt2(k),txres) end do draw(plot0) frame(wks) delete(res) wks = gsn_open_wks (output_type,"intensity_error_fcst_hour-" + g5nr_name(i) + "_" + start_date) res = True ; plot mods desired res@tiMainString = "Intensity Error (" + g5nr_name(i) + "): " + start_date res@trYMaxF = max_intensity + 5 res@trYMinF = 0 res@gsnFrame = False ; don't advance frame yet ;colors_plot = (/"Black","Red"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@vpWidthF = 0.8 ; set width and height res@vpHeightF = 0.3 res@vpXF = 0.15 res@vpYF = 0.9 res@trXMinF = 0 res@trXMaxF = x_axis_max x_axis = ispan(0,tointeger(x_axis_max),1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour res@tmXBValues = ispan(0,tointeger(x_axis_max),2) res@tmXBLabels = fcst_hr(0:tointeger(x_axis_max):2) res@tiYAxisString = "Intensity Error (hPa)" res@tiXAxisString = "Forecast lead time (hr)" res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.012 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 8.0 res@tmXMajorGrid = True ; implement x grid res@tmXMajorGridThicknessF = 1.0 ; 2.0 is default res@tmXMajorGridLineDashPattern = 2 res@tiMainFontHeightF = 0.021 res@xyLineColor = colors_plot(0) plot0 = gsn_csm_xy(wks,x_axis,intensity_error_ctl(i,check_ctl,0:x_axis_max),res) ; plot additional tracks res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 8.0 if (count_exp1(0).gt.count_ctl) then plot1 = gsn_csm_xy(wks,x_axis,intensity_error_exp(0,i,check_exp(0),0:x_axis_max),res) else plot1 = gsn_csm_xy(wks,x_axis(0:count_exp1(0)-1),intensity_error_exp(0,i,check_exp(0),0:count_exp1(0)-1),res) end if overlay(plot0,plot1) if (num_exp_cases.ge.2) then res@xyLineColor = colors_plot(2) res@xyLineThicknessF = 8.0 if (count_exp1(1).gt.count_ctl) then plot2 = gsn_csm_xy(wks,x_axis,intensity_error_exp(1,i,check_exp(1),0:x_axis_max),res) else plot2 = gsn_csm_xy(wks,x_axis(0:count_exp1(1)-1),intensity_error_exp(1,i,check_exp(1),0:count_exp1(1)-1),res) end if overlay(plot0,plot2) end if if (num_exp_cases.ge.3) then res@xyLineColor = colors_plot(3) res@xyLineThicknessF = 8.0 if (count_exp1(2).gt.count_ctl) then plot3 = gsn_csm_xy(wks,x_axis,intensity_error_exp(2,i,check_exp(2),0:x_axis_max),res) else plot3 = gsn_csm_xy(wks,x_axis(0:count_exp1(2)-1),intensity_error_exp(2,i,check_exp(2),0:count_exp1(2)-1),res) end if overlay(plot0,plot3) end if if (num_exp_cases.ge.4) then res@xyLineColor = colors_plot(4) res@xyLineThicknessF = 8.0 if (count_exp1(3).gt.count_ctl) then plot4 = gsn_csm_xy(wks,x_axis,intensity_error_exp(3,i,check_exp(3),0:x_axis_max),res) else plot4 = gsn_csm_xy(wks,x_axis(0:count_exp1(3)-1),intensity_error_exp(3,i,check_exp(3),0:count_exp1(3)-1),res) end if overlay(plot0,plot4) end if ; Attach a legend txres = True txres@txFontHeightF = 0.016 do k=0,dimsizes(labels2)-1 txres@txFontColor = colors_plot(k) gsn_text_ndc(wks,labels2(k),xtxt2(k),ytxt2(k),txres) end do draw(plot0) frame(wks) delete(res) delete(x_axis) ;;;; storm track map result ; determine bounds of track map maxes_lon = max(g5nr_lon(i,start_index:length_plot_g5)) maxes_lat = max(g5nr_lat(i,start_index:length_plot_g5)) mins_lon = min(g5nr_lon(i,start_index:length_plot_g5)) mins_lat = min(g5nr_lat(i,start_index:length_plot_g5)) ; used to get center longitude for pretty map spacing if (maxes_lon.lt.0) then maxes_lon = maxes_lon + 360.0 end if if (mins_lon.lt.0) then mins_lon = mins_lon + 360.0 end if wks=gsn_open_wks(output_type,"track_map_" + g5nr_name(i) + "_" + start_date) res = True res@gsnMaximize = True ; maximize plot in frame res@gsnFrame = False ; don't advance frame res@gsnDraw = False ; don't draw plot right away, wait till later res@mpFillOn = False res@mpPerimOn = True res@mpOutlineBoundarySets = "GeophysicalAndUSStates" res@pmTickMarkDisplayMode = "Always" res@mpLimitMode = "LatLon" ; select subregion res@mpMinLatF = mins_lat - 10.0 res@mpMaxLatF = maxes_lat + 10.0 res@mpMinLonF = mins_lon - 10.0 res@mpMaxLonF = maxes_lon + 10.0 res@mpCenterLonF = (mins_lon - 10.0) + (maxes_lon + 10.0) / 2.0 res@tmYROn = False ; turn off right and top tickmarks res@tmXTOn = False res@tiMainString = " " + g5nr_name(i) + ": " + start_date res@tiMainFontHeightF = 0.02 map = gsn_csm_map_ce(wks,res) ; create map ; Set up some legend resources. lgres = True lgres@lgLineColors = colors_plot(0:dimsizes(labels2)) lgres@lgLineThicknessF = 2. lgres@lgLabelFontHeightF = .08 ; set the legend label font thickness lgres@vpWidthF = 0.15 ; width of legend (NDC) lgres@vpHeightF = 0.1 ; height of legend (NDC) lgres@lgMonoDashIndex = True lgres@lgPerimColor = "orange" ; draw the box perimeter in orange lgres@lgPerimThicknessF = 2.0 ; thicken the box perimeter labels = array_append_record("G5NR",labels2,0) ; Create the legend. lbid = gsn_create_legend(wks,dimsizes(labels),labels,lgres) ; create legend ; Set up resources to attach legend to map. amres = True amres@amParallelPosF = -0.35 ;0.27 ; positive move legend to the right amres@amOrthogonalPosF = -0.35 ; positive move the legend down annoid1 = gsn_add_annotation(map,lbid,amres) ; attach legend to plot ; Add trajectory lines. pres = True ; polyline resources pres@gsLineThicknessF = 3.0 ; line thickness pres@gsLineColor = colors_plot(0) line1 = gsn_add_polyline(wks,map,g5nr_lon(i,start_index:length_plot_g5),g5nr_lat(i,start_index:length_plot_g5),pres) ; draw the traj pres = True ; polyline resources pres@gsLineColor = colors_plot(1) line2 = gsn_add_polyline(wks,map,ctl_lon_array(i,check_ctl,0:length_plot_ctl),ctl_lat_array(i,check_ctl,0:length_plot_ctl),pres) ; draw the traj pres = True ; polyline resources pres@gsLineColor = colors_plot(2) line3 = gsn_add_polyline(wks,map,exp_lon_array(0,i,check_exp(0),0:count_exp1(0)-1),exp_lat_array(0,i,check_exp(0),0:count_exp1(0)-1),pres) ; draw the traj if (num_exp_cases.ge.2) then pres = True ; polyline resources pres@gsLineColor = colors_plot(3) line4 = gsn_add_polyline(wks,map,exp_lon_array(1,i,check_exp(1),0:count_exp1(1)-1),exp_lat_array(1,i,check_exp(1),0:count_exp1(1)-1),pres) ; draw the traj end if if (num_exp_cases.ge.3) then pres = True ; polyline resources pres@gsLineColor = colors_plot(4) line5 = gsn_add_polyline(wks,map,exp_lon_array(2,i,check_exp(2),0:count_exp1(2)-1),exp_lat_array(2,i,check_exp(2),0:count_exp1(2)-1),pres) ; draw the traj end if if (num_exp_cases.ge.4) then pres = True ; polyline resources pres@gsLineColor = colors_plot(5) line6 = gsn_add_polyline(wks,map,exp_lon_array(3,i,check_exp(3),0:count_exp1(3)-1),exp_lat_array(3,i,check_exp(3),0:count_exp1(3)-1),pres) ; draw the traj end if ; Add markers to the trajectories. mres = True ; marker resources for best track mres@gsMarkerIndex = 16 ; marker style (filled circle) mres@gsMarkerSizeF = 8.0 ; marker size mres@gsMarkerColor = colors_plot(0) ; maker color markers1 = gsn_add_polymarker(wks,map,g5nr_lon(i,start_index:length_plot_g5),g5nr_lat(i,start_index:length_plot_g5),mres) mres@gsMarkerColor = colors_plot(1) ; maker color markers2 = gsn_add_polymarker(wks,map,ctl_lon_array(i,check_ctl,0:length_plot_ctl),ctl_lat_array(i,check_ctl,0:length_plot_ctl),mres) mres@gsMarkerColor = colors_plot(2) ; maker color markers3 = gsn_add_polymarker(wks,map,exp_lon_array(0,i,check_exp(0),0:count_exp1(0)-1),exp_lat_array(0,i,check_exp(0),0:count_exp1(0)-1),mres) if (num_exp_cases.ge.2) then mres@gsMarkerColor = colors_plot(3) ; maker color markers4 = gsn_add_polymarker(wks,map,exp_lon_array(1,i,check_exp(1),0:count_exp1(1)-1),exp_lat_array(1,i,check_exp(1),0:count_exp1(1)-1),mres) end if if (num_exp_cases.ge.3) then mres@gsMarkerColor = colors_plot(4) ; maker color markers5 = gsn_add_polymarker(wks,map,exp_lon_array(2,i,check_exp(2),0:count_exp1(2)-1),exp_lat_array(2,i,check_exp(2),0:count_exp1(2)-1),mres) end if if (num_exp_cases.ge.4) then mres@gsMarkerColor = colors_plot(5) ; maker color markers6 = gsn_add_polymarker(wks,map,exp_lon_array(3,i,check_exp(3),0:count_exp1(3)-1),exp_lat_array(3,i,check_exp(3),0:count_exp1(3)-1),mres) end if draw(map) frame(wks) delete(wks) delete(res) delete(mres) delete(pres) delete(lgres) delete(maxes_lon) delete(maxes_lat) delete(mins_lon) delete(mins_lat) delete(x_axis_max) delete(count_exp) delete(count_exp1) delete(x_axis_max_temp) end if end if ; end if for checking missing values end if end do ; end loop over GFS forecast runs delete(fcst_dates) delete(fcst_runs) else ; else for checking if both storms exist in NR and EXP print("Skipping storm: " + g5nr_name(i) + "...") end if ; end if for checking if storms exist in NR and EXP(s) end do ; end loop over computing stats over each storm of NR ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; average the track and intensity errors for both experiments print("Computing average track and intensity error statistics...") ;;;;;;;;;;;;;;;;;; loop over each control and experiment case and conform to each other so all have sample number of sample sizes ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; track_error_ctl_conform = new((/dimsizes(EXP_name),sizes_g5nr,40,29/),"float") wind_error_ctl_conform = new((/dimsizes(EXP_name),sizes_g5nr,40,29/),"float") intensity_error_ctl_conform = new((/dimsizes(EXP_name),sizes_g5nr,40,29/),"float") track_error_exp_conform = new((/dimsizes(EXP_name),sizes_g5nr,40,29/),"float") wind_error_exp_conform = new((/dimsizes(EXP_name),sizes_g5nr,40,29/),"float") intensity_error_exp_conform = new((/dimsizes(EXP_name),sizes_g5nr,40,29/),"float") do e=0,dimsizes(g5nr_name)-1 ; loop over G5NR storms do f=0,39 ; loop over GFS cycles n_ctl = num(.not.ismissing(track_error_ctl(e,f,:))) n_exp = new((dimsizes(EXP_name)),"integer") do a=0,dimsizes(EXP_name)-1 ; loop over experiment cases, conform each exp to ctl separately n_exp(a) = num(.not.ismissing(track_error_exp(a,e,f,:))) if (n_ctl.gt.0 .and. n_exp(a).gt.0) then ; means not all are missing for all cases appended1 = array_append_record(n_ctl,n_exp(a),0) ; combine to compute the min correct_length = min(appended1)-1 ; array indexing requires the minus one track_error_exp_conform(a,e,f,0:correct_length) = track_error_exp(a,e,f,0:correct_length) wind_error_exp_conform(a,e,f,0:correct_length) = wind_error_exp(a,e,f,0:correct_length) intensity_error_exp_conform(a,e,f,0:correct_length) = intensity_error_exp(a,e,f,0:correct_length) track_error_ctl_conform(a,e,f,0:correct_length) = track_error_ctl(e,f,0:correct_length) wind_error_ctl_conform(a,e,f,0:correct_length) = wind_error_ctl(e,f,0:correct_length) intensity_error_ctl_conform(a,e,f,0:correct_length) = intensity_error_ctl(e,f,0:correct_length) end if end do ; end loop over experiments end do ; end loop over GFS cycles end do ; end loop over G5NR storms ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; track_error_exp_conform!0 = "exp" track_error_exp_conform!1 = "storm" track_error_exp_conform!2 = "cycle" track_error_exp_conform!3 = "fcst_hr" intensity_error_exp_conform!0 = "exp" intensity_error_exp_conform!1 = "storm" intensity_error_exp_conform!2 = "cycle" intensity_error_exp_conform!3 = "fcst_hr" wind_error_exp_conform!0 = "exp" wind_error_exp_conform!1 = "storm" wind_error_exp_conform!2 = "cycle" wind_error_exp_conform!3 = "fcst_hr" track_error_ctl_conform!0 = "exp" track_error_ctl_conform!1 = "storm" track_error_ctl_conform!2 = "cycle" track_error_ctl_conform!3 = "fcst_hr" intensity_error_ctl_conform!0 = "exp" intensity_error_ctl_conform!1 = "storm" intensity_error_ctl_conform!2 = "cycle" intensity_error_ctl_conform!3 = "fcst_hr" wind_error_ctl_conform!0 = "exp" wind_error_ctl_conform!1 = "storm" wind_error_ctl_conform!2 = "cycle" wind_error_ctl_conform!3 = "fcst_hr" track_error_exp_reorder = track_error_exp_conform(storm|:,fcst_hr|:,exp|:,cycle|:) intensity_error_exp_reorder = intensity_error_exp_conform(storm|:,fcst_hr|:,exp|:,cycle|:) wind_error_exp_reorder = wind_error_exp_conform(storm|:,fcst_hr|:,exp|:,cycle|:) track_error_ctl_reorder = track_error_ctl_conform(storm|:,fcst_hr|:,exp|:,cycle|:) intensity_error_ctl_reorder = intensity_error_ctl_conform(storm|:,fcst_hr|:,exp|:,cycle|:) wind_error_ctl_reorder = wind_error_ctl_conform(storm|:,fcst_hr|:,exp|:,cycle|:) track_error_ctl_st_fh = dim_avg_n(track_error_ctl_reorder,(/3/)) ; average as function of storm and forecast hour track_error_exp_st_fh = dim_avg_n(track_error_exp_reorder,(/3/)) intensity_error_ctl_st_fh = dim_avg_n(intensity_error_ctl_reorder,(/3/)) intensity_error_exp_st_fh = dim_avg_n(intensity_error_exp_reorder,(/3/)) wind_error_ctl_st_fh = dim_avg_n(wind_error_ctl_reorder,(/3/)) wind_error_exp_st_fh = dim_avg_n(wind_error_exp_reorder,(/3/)) track_error_ctl_fh = dim_avg_n(track_error_ctl_conform,(/1,2/)) ; average as function of forecast hour only track_error_exp_fh = dim_avg_n(track_error_exp_conform,(/1,2/)) track_error_ctl_fh_var = dim_variance_n(track_error_ctl_conform,(/1,2/)) ; average as function of forecast hour only track_error_exp_fh_var = dim_variance_n(track_error_exp_conform,(/1,2/)) track_error_ctl_fh_std = dim_stddev_n(track_error_ctl_conform,(/1,2/)) ; average as function of forecast hour only track_error_exp_fh_std = dim_stddev_n(track_error_exp_conform,(/1,2/)) intensity_error_ctl_fh_var = dim_variance_n(intensity_error_ctl_conform,(/1,2/)) ; average as function of forecast hour only intensity_error_exp_fh_var = dim_variance_n(intensity_error_exp_conform,(/1,2/)) intensity_error_ctl_fh_std = dim_stddev_n(intensity_error_ctl_conform,(/1,2/)) ; average as function of forecast hour only intensity_error_exp_fh_std = dim_stddev_n(intensity_error_exp_conform,(/1,2/)) intensity_error_ctl_fh = dim_avg_n(intensity_error_ctl_conform,(/1,2/)) intensity_error_exp_fh = dim_avg_n(intensity_error_exp_conform,(/1,2/)) wind_error_ctl_fh_var = dim_variance_n(wind_error_ctl_conform,(/1,2/)) ; average as function of forecast hour only wind_error_exp_fh_var = dim_variance_n(wind_error_exp_conform,(/1,2/)) wind_error_ctl_fh_std = dim_stddev_n(wind_error_ctl_conform,(/1,2/)) ; average as function of forecast hour only wind_error_exp_fh_std = dim_stddev_n(wind_error_exp_conform,(/1,2/)) wind_error_ctl_fh = dim_avg_n(wind_error_ctl_conform,(/1,2/)) wind_error_exp_fh = dim_avg_n(wind_error_exp_conform,(/1,2/)) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; find the longest length ctl based on its comparison with each experiment so that we plot the longest ctl max_ctl_length = dim_num_n(.not.ismissing(track_error_ctl_fh),1) ; max length for each paired experiment location_max_temp = ind(max_ctl_length.eq.max(max_ctl_length)) ; there could be more than one index due to matching lengths, so we always choose the first index when plotting location_max = location_max_temp(0) ; total average over all storms and hours, one value here track_error_ctl_avg = avg(track_error_ctl_conform) ; total average track_error_exp_avg = dim_avg_n(track_error_exp_conform,(/1,2,3/)) intensity_error_ctl_avg = avg(intensity_error_ctl_conform) intensity_error_exp_avg = dim_avg_n(intensity_error_exp_conform,(/1,2,3/)) wind_error_ctl_avg = avg(wind_error_ctl_conform) wind_error_exp_avg = dim_avg_n(wind_error_exp_conform,(/1,2,3/)) ; determine number of non-missing values at each forecast hour in each experiment do a=0,dimsizes(EXP_name)-1 do i=0,28 ; loop over forecast hours (output every 6 hours) fcst_hr_array_exp(a,i) = num(.not.ismissing(track_error_exp_conform(a,:,:,i))) fcst_hr_array_ctl(a,i) = num(.not.ismissing(track_error_ctl_conform(a,:,:,i))) end do end do ; compute significance over forecast hour for track and intensity and wind ; loop over the experiments and compare to CTL prob_fh_track = new((/dimsizes(EXP_name),29/),"float") prob_fh_intensity = prob_fh_track prob_fh_wind = prob_fh_track track_error_diff_std = prob_fh_track intensity_error_diff_std = track_error_diff_std wind_error_diff_std = track_error_diff_std interval_track_diff = prob_fh_track interval_intensity_diff = prob_fh_track interval_wind_diff = prob_fh_track track_error_diff_avg = prob_fh_track intensity_error_diff_avg = prob_fh_track wind_error_diff_avg = prob_fh_track number_sample = prob_fh_track const_z = new((/29/),"float") ; 95% significance level for paired t test do a=0,dimsizes(EXP_name)-1 track_error_diff_std(a,:) = dim_stddev_n(track_error_exp_conform(a,:,:,:)-track_error_ctl_conform(a,:,:,:),(/0,1/)) ; forecast hour only intensity_error_diff_std(a,:) = dim_stddev_n(intensity_error_exp_conform(a,:,:,:)-intensity_error_ctl_conform(a,:,:,:),(/0,1/)) ; forecast hour only, wind_error_diff_std(a,:) = dim_stddev_n(wind_error_exp_conform(a,:,:,:)-wind_error_ctl_conform(a,:,:,:),(/0,1/)) ; forecast hour only track_error_diff_avg(a,:) = dim_avg_n(track_error_exp_conform(a,:,:,:)-track_error_ctl_conform(a,:,:,:),(/0,1/)) ; forecast hour only intensity_error_diff_avg(a,:) = dim_avg_n(intensity_error_exp_conform(a,:,:,:)-intensity_error_ctl_conform(a,:,:,:),(/0,1/)) ; forecast hour only wind_error_diff_avg(a,:) = dim_avg_n(wind_error_exp_conform(a,:,:,:)-wind_error_ctl_conform(a,:,:,:),(/0,1/)) ; forecast hour only, do y=0,28 ; loop over samples to divide by smallest of the two samples if (fcst_hr_array_ctl(a,y).le.fcst_hr_array_exp(a,y)) then number_sample(a,y) = fcst_hr_array_ctl(a,y) if (number_sample(a,y).lt.20) then const_z(y) = 2.228 end if if (number_sample(a,y).ge.20 .and. number_sample(a,y).lt.40) then const_z(y) = 2.042 end if if (number_sample(a,y).ge.40 .and. number_sample(a,y).lt.80) then const_z(y) = 2.000 end if if (number_sample(a,y).ge.80) then const_z(y) = 1.96 end if else number_sample(a,y) = fcst_hr_array_exp(a,y) if (number_sample(a,y).lt.20) then const_z(y) = 2.228 end if if (number_sample(a,y).ge.20 .and. number_sample(a,y).lt.40) then const_z(y) = 2.042 end if if (number_sample(a,y).ge.40 .and. number_sample(a,y).lt.80) then const_z(y) = 2.000 end if if (number_sample(a,y).ge.80) then const_z(y) = 1.96 end if end if end do number_valid = ind(number_sample(a,:).ne.0) if (dimsizes(number_valid).ge.1) then interval_track_diff(a,number_valid) = const_z(number_valid)*track_error_diff_std(a,number_valid) / sqrt(number_sample(a,number_valid)) ; % 95 interval interval_intensity_diff(a,number_valid) = const_z(number_valid)*intensity_error_diff_std(a,number_valid) / sqrt(number_sample(a,number_valid)) ; 95% interval interval_wind_diff(a,number_valid) = const_z(number_valid)*wind_error_diff_std(a,number_valid) / sqrt(number_sample(a,number_valid)) ; 95% interval end if delete(number_valid) end do ; end loop over EXP names ; plot results print("Plotting results...") ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Storm track errors as histogram, along with wind and intensity errors ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;colors_plot = (/"Black","Red","Orange","Green","Blue","Yellow"/) if (num_exp_cases.eq.1) then names_short_1 = (/"A","B","","A","B"/) ; names for legend names_short_2 = (/"A = "+CTL_name,"B = "+EXP_name(0)/) ; names for legend names_legend = (/"A","B"/) end if if (num_exp_cases.eq.2) then names_short_1 = (/"A","B","C","","A","B","C"/) ; names for legend names_short_2 = (/"A = "+CTL_name,"B = "+EXP_name(0),"C = "+EXP_name(1)/) ; names for legend names_legend = (/"A","B","C"/) end if if (num_exp_cases.eq.3) then names_short_1 = (/"A","B","C","D","","A","B","C","D"/) ; names for legend names_short_2 = (/"A = "+CTL_name,"B = "+EXP_name(0),"C = "+EXP_name(1),"D = "+EXP_name(2)/) ; names for legend names_legend = (/"A","B","C","D"/) end if if (num_exp_cases.eq.4) then names_short_1 = (/"A","B","C","D","E","","A","B","C","D","E"/) ; names for legend names_short_2 = (/"A = "+CTL_name,"B = "+EXP_name(0),"C = "+EXP_name(1),"D = "+EXP_name(2),"E = "+EXP_name(3)/) ; names for legend names_legend = (/"A","B","C","D","E"/) end if xtxt = (/0.25,0.25,0.25,0.25,0.25,0.25,0.25,0.25/) ytxt = (/0.5,0.47,0.44,0.41,0.38,0.35,0.32,0.29/) print("Making track, wind, and intensity error plots as histograms...") if (num_exp_cases.eq.1) then x_axis_1 = ispan(2,4,2) x_axis_2 = ispan(8,10,2) y_axis_values_1 = (/track_error_ctl_avg,track_error_exp_avg(0)/) ; y values for the bar charts y_axis_values_2 = (/intensity_error_ctl_avg,intensity_error_exp_avg(0)/) y_axis_values_3 = (/wind_error_ctl_avg,wind_error_exp_avg(0)/) ; y values for the bar charts end if if (num_exp_cases.eq.2) then x_axis_1 = ispan(2,6,2) x_axis_2 = ispan(10,14,2) y_axis_values_1 = (/track_error_ctl_avg,track_error_exp_avg(0),track_error_exp_avg(1)/) ; y values for the bar charts y_axis_values_2 = (/intensity_error_ctl_avg,intensity_error_exp_avg(0),intensity_error_exp_avg(1)/) y_axis_values_3 = (/wind_error_ctl_avg,wind_error_exp_avg(0),wind_error_exp_avg(1)/) ; y values for the bar charts end if if (num_exp_cases.eq.3) then x_axis_1 = ispan(2,8,2) x_axis_2 = ispan(12,18,2) y_axis_values_1 = (/track_error_ctl_avg,track_error_exp_avg(0),track_error_exp_avg(1),track_error_exp_avg(2)/) ; y values for the bar charts y_axis_values_2 = (/intensity_error_ctl_avg,intensity_error_exp_avg(0),intensity_error_exp_avg(1),intensity_error_exp_avg(2)/) y_axis_values_3 = (/wind_error_ctl_avg,wind_error_exp_avg(0),wind_error_exp_avg(1),wind_error_exp_avg(2)/) ; y values for the bar charts end if if (num_exp_cases.eq.4) then x_axis_1 = ispan(2,10,2) x_axis_2 = ispan(14,22,2) y_axis_values_1 = (/track_error_ctl_avg,track_error_exp_avg(0),track_error_exp_avg(1),track_error_exp_avg(2),track_error_exp_avg(3)/) ; y values for the bar charts y_axis_values_2 = (/intensity_error_ctl_avg,intensity_error_exp_avg(0),intensity_error_exp_avg(1),intensity_error_exp_avg(2),intensity_error_exp_avg(3)/) y_axis_values_3 = (/wind_error_ctl_avg,wind_error_exp_avg(0),wind_error_exp_avg(1),wind_error_exp_avg(2),wind_error_exp_avg(3)/) ; y values for the bar charts end if max_y1=max(y_axis_values_1) max_y2=max(y_axis_values_2) min_y1=min(y_axis_values_1) min_y2=min(y_axis_values_2) wks = gsn_open_wks(output_type,"track-intensity-error-all-storms_" + region_file) ; create workstation res = True res@gsnFrame = False res@gsnDraw = False res@tiMainString = "Track and Intensity Error (" + count_storms + " TCs)" + " " + region_name res@gsnXYBarChart = True ; Create bar plot ;res@tiMainOffsetYF = 0.006 res@tiMainFontHeightF = 0.02 res@tmXTMode = "Explicit" ; Define own tick mark labels. res@tmXTValues = x_axis_1 res@tmXTOn = False ; turn on the top tickmarks res@tmXBOn = True res@tmXBLabelsOn = True res@tmXBLabelJust = "BottomCenter" res@gsnScale = True res@tmYROn = False res@tmYRMinorOn = True res@tmYLOn = True res@tmYLMinorOn = True res@gsnXYBarChartOutlineThicknessF = 1.5 res@tmXTMinorOn = False ;res@tiYAxisOffsetXF = 0.035 res@tiYAxisString = "Track Error (km)" res@tiYAxisFontHeightF = 0.018 res@trYMaxF = max_y1 + 200.0 ; max value on y-axis res@trYMinF = 0 res@gsnXYBarChartBarWidth = 1.5 ; change bar widths 0.75 originally res@gsnXYBarChartColors = colors_plot(1) res@tmXBMode = "Explicit" res@trXMinF = 0 ; adds space on either end res@trXMaxF = max(x_axis_2)+2.0 ; of the 1st and last bars res@tmXBValues = x_axis_1 res@tmXBLabels = names_legend res@tmXBLabelJust = "BottomCenter" plot = gsn_csm_xy(wks,x_axis_1,y_axis_values_1,res) bres = True bres@tmXBMode = "Explicit" bres@gsnFrame = False bres@gsnDraw = False bres@gsnXYBarChart = True bres@gsnXYBarChartColors = colors_plot(2) bres@tmYROn = True bres@tmYRMinorOn = True bres@tmYLOn = False bres@trYMinF = min_y2 - 5.0 ; hPa bres@tmXBValues = x_axis_2 bres@tiMainString = "" bres@trYMaxF = max_y2 + 5.0 ; hPa bres@tmYRLabelsOn = True bres@tiYAxisFontHeightF = 0.018 bres@gsnXYBarChartOutlineThicknessF = 1.5 bres@tmYROn = True bres@gsnXYBarChartBarWidth = 1.5 ; change bar widths 0.75 originally bres@tmXBOn = True bres@trXMinF = 0 ; adds space on either end bres@trXMaxF = max(x_axis_2)+2.0 ; of the 1st and last bars bres@tiYAxisString = "Intensity Error (hPa)" bres@tiYAxisSide = "Right" bres@tiYAxisAngleF = 270 bres@tmXBLabels = names_legend bres@tmXBLabelJust = "BottomCenter" bres@tmXTOn = False draw(plot) ;********************************************************** ; add text labels ;********************************************************** txres = True ; text mods desired txres@txFontHeightF = 0.016 ; default size is HUGE! txres@txAngleF = 0 ; text angle txres@txJust = "CenterCenter" ; puts text on top of bars; txres@txFontColor = "black" if (num_exp_cases.eq.1) then gsn_text_ndc(wks,names_short_2(0),0.35,0.75,txres) ; add labels gsn_text_ndc(wks,names_short_2(1),0.35,0.72,txres) ; add labels end if if (num_exp_cases.eq.2) then gsn_text_ndc(wks,names_short_2(0),0.35,0.75,txres) ; add labels gsn_text_ndc(wks,names_short_2(1),0.35,0.72,txres) ; add labels gsn_text_ndc(wks,names_short_2(2),0.35,0.69,txres) ; add labels end if if (num_exp_cases.eq.3) then gsn_text_ndc(wks,names_short_2(0),0.35,0.75,txres) ; add labels gsn_text_ndc(wks,names_short_2(1),0.35,0.72,txres) ; add labels gsn_text_ndc(wks,names_short_2(2),0.35,0.69,txres) ; add labels gsn_text_ndc(wks,names_short_2(3),0.35,0.66,txres) ; add labels end if if (num_exp_cases.eq.4) then gsn_text_ndc(wks,names_short_2(0),0.35,0.75,txres) ; add labels gsn_text_ndc(wks,names_short_2(1),0.35,0.72,txres) ; add labels gsn_text_ndc(wks,names_short_2(2),0.35,0.69,txres) ; add labels gsn_text_ndc(wks,names_short_2(3),0.35,0.66,txres) ; add labels gsn_text_ndc(wks,names_short_2(4),0.35,0.63,txres) ; add labels end if plot2 = gsn_csm_xy(wks,x_axis_2,y_axis_values_2,bres) ;---Create a blank plot in order to get additional X tickmark labels vres = True vres@tmXBMode = "Explicit" vres@tmXBValues = (/avg(x_axis_1),avg(x_axis_2)/) vres@tmXBLabels = (/"Position Error","Magnitude Error"/) vres@tmXBLabelDeltaF = 0.5 vres@tmXBLabelFontHeightF = 0.018 vres@tmYLOn = False vres@tmYROn = False ; vres@tmXBOn = False vres@tmXTOn = False blank_plot = gsn_csm_blank_plot(wks,vres) overlay(plot2,blank_plot) draw(plot2) frame(wks) delete(bres) delete(res) delete(vres) delete(txres) ; histogram of wind, keep same settings as prior histogram to avoid rewritting more code max_y1=max(y_axis_values_3) ; copied from above min_y1=min(y_axis_values_3) wks = gsn_open_wks(output_type,"wind-intensity-error-all-storms_" + region_file) ; create workstation res = True res@gsnFrame = False res@gsnDraw = False res@tiMainString = "Wind Error (" + count_storms + " TCs)" + " " + region_name res@gsnXYBarChart = True ; Create bar plot ;res@tiMainOffsetYF = 0.006 res@tiMainFontHeightF = 0.02 res@tmXTMode = "Explicit" ; Define own tick mark labels. res@tmXTValues = x_axis_1 res@tmXTOn = False ; turn on the top tickmarks res@tmXBOn = True res@tmXBLabelsOn = True res@tmXBLabelJust = "BottomCenter" res@gsnScale = True res@tmYROn = False res@tmYRMinorOn = True res@tmYLOn = True res@tmYLMinorOn = True res@gsnXYBarChartOutlineThicknessF = 1.5 res@tmXTMinorOn = False ;res@tiYAxisOffsetXF = 0.035 res@tiYAxisString = "Wind Error (kts)" res@tiYAxisFontHeightF = 0.018 res@trYMaxF = max_y1 + 15.0 ; max value on y-axis res@trYMinF = 0 res@gsnXYBarChartBarWidth = 1.5 ; change bar widths 0.75 originally res@gsnXYBarChartColors = colors_plot(1) res@tmXBMode = "Explicit" res@trXMinF = 0 ; adds space on either end res@trXMaxF = max(x_axis_1)+2.0 ; of the 1st and last bars res@tmXBValues = x_axis_1 res@tmXBLabels = names_legend res@tmXBLabelJust = "BottomCenter" plot = gsn_csm_xy(wks,x_axis_1,y_axis_values_3,res) draw(plot) ;********************************************************** ; add text labels ;********************************************************** txres = True ; text mods desired txres@txFontHeightF = 0.016 ; default size is HUGE! txres@txAngleF = 0 ; text angle txres@txJust = "CenterCenter" ; puts text on top of bars; txres@txFontColor = "black" if (num_exp_cases.eq.1) then gsn_text_ndc(wks,names_short_2(0),0.35,0.75,txres) ; add labels gsn_text_ndc(wks,names_short_2(1),0.35,0.72,txres) ; add labels end if if (num_exp_cases.eq.2) then gsn_text_ndc(wks,names_short_2(0),0.35,0.75,txres) ; add labels gsn_text_ndc(wks,names_short_2(1),0.35,0.72,txres) ; add labels gsn_text_ndc(wks,names_short_2(2),0.35,0.69,txres) ; add labels end if if (num_exp_cases.eq.3) then gsn_text_ndc(wks,names_short_2(0),0.35,0.75,txres) ; add labels gsn_text_ndc(wks,names_short_2(1),0.35,0.72,txres) ; add labels gsn_text_ndc(wks,names_short_2(2),0.35,0.69,txres) ; add labels gsn_text_ndc(wks,names_short_2(3),0.35,0.66,txres) ; add labels end if if (num_exp_cases.eq.4) then gsn_text_ndc(wks,names_short_2(0),0.35,0.75,txres) ; add labels gsn_text_ndc(wks,names_short_2(1),0.35,0.72,txres) ; add labels gsn_text_ndc(wks,names_short_2(2),0.35,0.69,txres) ; add labels gsn_text_ndc(wks,names_short_2(3),0.35,0.66,txres) ; add labels gsn_text_ndc(wks,names_short_2(4),0.35,0.63,txres) ; add labels end if frame(wks) ;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; plot errors as a function of forecast lead time ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; delete(res) delete(txres) ; determine max and min for track and intensity values ; pad both sides track_temp = array_append_record(track_error_ctl_fh(location_max,:),ndtooned(track_error_exp_fh),0) max_track = max(track_temp) min_track = min(track_temp) int_temp = array_append_record(intensity_error_ctl_fh(location_max,:),ndtooned(intensity_error_exp_fh),0) max_int = max(int_temp) min_int = min(int_temp) wind_temp = array_append_record(wind_error_ctl_fh(location_max,:),ndtooned(wind_error_exp_fh),0) max_wind = max(wind_temp) min_wind = min(wind_temp) max_fcst = max(fcst_hr_array_ctl) min_fcst = min(fcst_hr_array_ctl) wks = gsn_open_wks (output_type,"track_error_fcst_hour-all-storms_" + region_file) res = True ; plot mods desired res@gsnXYBarChart = False res@tiMainString = "Track Error (" + count_storms + " TCs)" + " " + region_name res@trYMaxF = max_track + 200.0 res@trYMinF = 0 res@gsnFrame = False ; don't advance frame yet res@gsnDraw = False res@gsnFrame = False res@vpXF = 0.15 res@vpYF = 0.90 res@vpWidthF = 0.7 res@vpHeightF = 0.6 res@gsnScale = True res@trXMinF = 0 res@trXMaxF = 28 x_axis = ispan(0,28,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour res@tmXBValues = ispan(0,28,2) res@tmXBLabels = ispan(0,168,12) res@tiYAxisString = "Track Error (km)" res@tiXAxisString = "Forecast lead time (hr)" res@tiXAxisFontHeightF = 0.015 res@tmYROn = False res@tmYROn = False res@tmYRMinorOn = False res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 8.0 res@tmXMajorGrid = True ; implement x grid res@tmXMajorGridThicknessF = 1.0 ; 2.0 is default res@tmXMajorGridLineDashPattern = 2 res@tiMainFontHeightF = 0.021 res@xyLineColor = colors_plot(0) plot0 = gsn_csm_xy(wks,x_axis,track_error_ctl_fh(location_max,0:28),res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 8.0 plot1 = gsn_csm_xy(wks,x_axis,track_error_exp_fh(0,0:28),res) if (num_exp_cases.ge.2) then res@xyLineColor = colors_plot(2) res@xyLineThicknessF = 8.0 plot2 = gsn_csm_xy(wks,x_axis,track_error_exp_fh(1,0:28),res) overlay(plot0,plot2) end if if (num_exp_cases.ge.3) then res@xyLineColor = colors_plot(3) res@xyLineThicknessF = 8.0 plot3 = gsn_csm_xy(wks,x_axis,track_error_exp_fh(2,0:28),res) overlay(plot0,plot3) end if if (num_exp_cases.ge.4) then res@xyLineColor = colors_plot(4) res@xyLineThicknessF = 8.0 plot4 = gsn_csm_xy(wks,x_axis,track_error_exp_fh(3,0:28),res) overlay(plot0,plot4) end if bres = True bres@tmXBMode = "Explicit" bres@gsnFrame = False bres@gsnDraw = False bres@vpXF = 0.15 bres@vpYF = 0.90 bres@vpWidthF = 0.7 bres@vpHeightF = 0.6 bres@gsnXYBarChart = True bres@gsnXYBarChartColors = "grey78" bres@tmYROn = True bres@tmYRMinorOn = True bres@tmYLOn = False bres@trYMinF = 0 bres@tiMainString = "" bres@trYMaxF = max_fcst + 4.0 ; number of samples in each forecast hour bres@tmYRLabelsOn = True bres@tiYAxisFontHeightF = 0.015 bres@tmYRLabelFontHeightF = 0.015 bres@gsnXYBarChartOutlineThicknessF = 1.5 bres@tmYROn = True bres@gsnXYBarChartBarWidth = 0.5 ; change bar widths 0.75 originally bres@tmXBOn = False bres@trXMinF = 0 ; adds space on either end bres@trXMaxF = 28 ; of the 1st and last bars bres@tiYAxisString = "# fcsts" bres@tiYAxisSide = "Right" bres@tiYAxisAngleF = 270 bres@tmXTOn = False plot5 = gsn_csm_xy(wks,x_axis,fcst_hr_array_ctl(location_max,0:28),bres) bres@gsnXYBarChartColors = "grey78" ;grey71 ;plot6 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp(0,0:28),bres) ;draw(plot5) ;draw(plot6) if (num_exp_cases.ge.2) then bres@gsnXYBarChartColors = "grey78" ; grey65 ;plot7 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp(1,0:28),bres) ;draw(plot7) end if if (num_exp_cases.ge.3) then bres@gsnXYBarChartColors = "grey78" ; grey59 ;plot8 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp(2,0:28),bres) ;draw(plot8) end if if (num_exp_cases.ge.4) then bres@gsnXYBarChartColors = "grey78" ; grey58 ;plot9 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp(3,0:28),bres) ;draw(plot9) end if overlay(plot0,plot1) ; Attach a legend marker_value = (/3,4,6,9/) mres = True txres = True txres@txFontHeightF = 0.016 do k=0,dimsizes(labels2)-1 txres@txFontColor = colors_plot(k) gsn_text_ndc(wks,labels2(k),xtxt3(k),ytxt3(k),txres) if (k.ge.1) then mres@gsMarkerIndex = marker_value(k-1) mres@gsMarkerSizeF = 16.0 mres@gsMarkerThicknessF = 4.0 mres@gsMarkerColor = colors_plot(k) gsn_polymarker_ndc(wks,xtxt3(k)-0.07,ytxt3(k),mres) end if end do ; attach significance do k=0,dimsizes(EXP_name)-1 do w=0,28 ; loop over forecast hours mres@gsMarkerIndex = marker_value(k) mres@gsMarkerSizeF = 16.0 mres@gsMarkerThicknessF = 4.0 ;mres@gsMarkerColor = colors_plot(k+1) if (.not.ismissing(track_error_diff_avg(k,w)) .and. interval_track_diff(k,w).gt.0) then if (abs(track_error_diff_avg(k,w)).gt.interval_track_diff(k,w)) then ; 95% confidence interval if (track_error_diff_avg(k,w).lt.0) then ; then exp is better mres@gsMarkerColor = colors_plot(k+1) else ; ctl is better mres@gsMarkerColor = colors_plot(0) end if gsn_polymarker(wks,plot0,x_axis(w),20.0,mres) end if end if end do end do ;;;;; resbot = True resbot@vpXF = 0.15 ; The left side of the box location resbot@vpYF = 0.20 ; The top side of the plot box location resbot@vpWidthF = 0.70 ; The Width of the plot box resbot@vpHeightF = 0.05 ; The height of the plot box resbot@trYMaxF = 3 ; random value does not need to be changed resbot@trXMinF = 0 resbot@trXMaxF = 28 resbot@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour resbot@tmXBValues = ispan(0,28,2) resbot@tmXBLabels = ispan(0,168,12) resbot@tiYAxisString = "# fcsts" resbot@tmYLBorderOn = False resbot@tmYLMinorPerMajor = 1 resbot@tmYLOn = False resbot@tmYRBorderOn = False resbot@tmYROn = False resbot@tiXAxisString = "Forecast lead time (hr)" resbot@tiXAxisFontHeightF = 0.015 resbot@tmYROn = False resbot@tmYROn = False resbot@tmYRMinorOn = False resbot@tiYAxisFontHeightF = 0.015 resbot@tmXBLabelFontHeightF = 0.015 resbot@tmYLLabelFontHeightF = 0.015 resbot@xyLineThicknessF = 0.05 resbot@tiMainFontHeightF = 0.021 resbot@tmXTOn = False ; turn off the top tickmarks resbot@tmYLOn = False resbot@tmXBMinorOn = False ; No minor tick marks. resbot@gsnDraw = False ; Don't draw individual plot. resbot@gsnFrame = False ; Don't advance frame. ;resbot@gsnYRefLine = 0.0 ; create a reference line ;resbot@gsnAboveYRefLineColor = "green" ; above ref line fill green resbot@gsnXYBarChart = False ; turn on bar chart plot_blank = gsn_csm_xy(wks,x_axis,0.0*fcst_hr_array_ctl(location_max,0:28),resbot) ; plot false values to create plot, then overlay text values on plot ; attache the number of paired cases over the plot y_value_text = (/1.5,1.5,1.5/) do k=0,dimsizes(EXP_name)-1 do w=0,28,2 ; loop over forecast hours, skip first hour for simplicity of viewing image qres = True qres@txFontHeightF = 0.017 qres@txJust = "CenterCenter" qres@txFontColor = colors_plot(k+1) ; color coded by experiment text = gsn_add_text(wks,plot_blank,str_squeeze(tostring(tointeger(fcst_hr_array_exp(k,w)))),x_axis(w),y_value_text(k),qres) end do end do ;;;;;;;;;;;;; draw(plot0) draw(plot_blank) frame(wks) delete(res) delete(bres) delete(resbot) ; wind wks = gsn_open_wks (output_type,"wind_error_fcst_hour-all-storms_" + region_file) res = True ; plot mods desired res@gsnXYBarChart = False res@tiMainString = "Wind Error (" + count_storms + " TCs)" + " " + region_name res@trYMaxF = max_wind + 20.0 res@trYMinF = 0 res@gsnFrame = False ; don't advance frame yet res@gsnDraw = False res@gsnFrame = False res@vpXF = 0.15 res@vpYF = 0.90 res@vpWidthF = 0.7 res@vpHeightF = 0.6 res@gsnScale = True res@trXMinF = 0 res@trXMaxF = 28 x_axis = ispan(0,28,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour res@tmXBValues = ispan(0,28,2) res@tmXBLabels = ispan(0,168,12) res@tiYAxisString = "Wind Error (kts)" res@tiXAxisString = "Forecast lead time (hr)" res@tiXAxisFontHeightF = 0.015 res@tmYROn = False res@tmYROn = False res@tmYRMinorOn = False res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@tmXMajorGrid = True ; implement x grid res@tmXMajorGridThicknessF = 1.0 ; 2.0 is default res@tmXMajorGridLineDashPattern = 2 res@xyLineThicknessF = 8.0 res@tiMainFontHeightF = 0.021 res@xyLineColor = colors_plot(0) plot0 = gsn_csm_xy(wks,x_axis,wind_error_ctl_fh(location_max,0:28),res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 8.0 plot1 = gsn_csm_xy(wks,x_axis,wind_error_exp_fh(0,0:28),res) if (num_exp_cases.ge.2) then res@xyLineColor = colors_plot(2) res@xyLineThicknessF = 8.0 plot2 = gsn_csm_xy(wks,x_axis,wind_error_exp_fh(1,0:28),res) overlay(plot0,plot2) end if if (num_exp_cases.ge.3) then res@xyLineColor = colors_plot(3) res@xyLineThicknessF = 8.0 plot3 = gsn_csm_xy(wks,x_axis,wind_error_exp_fh(2,0:28),res) overlay(plot0,plot3) end if if (num_exp_cases.ge.4) then res@xyLineColor = colors_plot(4) res@xyLineThicknessF = 8.0 plot4 = gsn_csm_xy(wks,x_axis,wind_error_exp_fh(3,0:28),res) overlay(plot0,plot4) end if bres = True bres@tmXBMode = "Explicit" bres@gsnFrame = False bres@gsnDraw = False bres@gsnXYBarChart = True bres@gsnXYBarChartColors = "grey78" bres@tmYROn = True bres@tmYRMinorOn = True bres@vpXF = 0.15 bres@vpYF = 0.90 bres@vpWidthF = 0.7 bres@vpHeightF = 0.6 bres@tmYLOn = False bres@trYMinF = 0 bres@tiMainString = "" bres@trYMaxF = max_fcst + 4.0 ; number of samples in each forecast hour bres@tmYRLabelsOn = True bres@tiYAxisFontHeightF = 0.015 bres@tmYRLabelFontHeightF = 0.015 bres@gsnXYBarChartOutlineThicknessF = 1.5 bres@tmYROn = True bres@gsnXYBarChartBarWidth = 0.5 ; change bar widths 0.75 originally bres@tmXBOn = False bres@trXMinF = 0 ; adds space on either end bres@trXMaxF = 28 ; of the 1st and last bars bres@tiYAxisString = "# fcsts" bres@tiYAxisSide = "Right" bres@tiYAxisAngleF = 270 bres@tmXTOn = False plot5 = gsn_csm_xy(wks,x_axis,fcst_hr_array_ctl(location_max,0:28),bres) bres@gsnXYBarChartColors = "grey78" ;plot6 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp(0,0:28),bres) ;draw(plot5) ;draw(plot6) if (num_exp_cases.ge.2) then bres@gsnXYBarChartColors = "grey78" ;plot7 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp(1,0:28),bres) ;draw(plot7) end if if (num_exp_cases.ge.3) then bres@gsnXYBarChartColors = "grey78" ;plot8 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp(2,0:28),bres) ;draw(plot8) end if if (num_exp_cases.ge.4) then bres@gsnXYBarChartColors = "grey78" ;plot9 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp(3,0:28),bres) ;draw(plot9) end if overlay(plot0,plot1) ; Attach a legend marker_value = (/3,4,6,9/) mres = True txres = True txres@txFontHeightF = 0.016 do k=0,dimsizes(labels2)-1 txres@txFontColor = colors_plot(k) gsn_text_ndc(wks,labels2(k),xtxt3(k),ytxt3(k),txres) if (k.ge.1) then mres@gsMarkerIndex = marker_value(k-1) mres@gsMarkerSizeF = 16.0 mres@gsMarkerThicknessF = 4.0 mres@gsMarkerColor = colors_plot(k) gsn_polymarker_ndc(wks,xtxt3(k)-0.07,ytxt3(k),mres) end if end do ; attach significance do k=0,dimsizes(EXP_name)-1 do w=0,28 ; loop over forecast hours mres@gsMarkerIndex = marker_value(k) mres@gsMarkerSizeF = 16.0 mres@gsMarkerThicknessF = 4.0 ;mres@gsMarkerColor = colors_plot(k+1) if (.not.ismissing(wind_error_diff_avg(k,w)) .and. interval_wind_diff(k,w).gt.0) then if (abs(wind_error_diff_avg(k,w)).gt.interval_wind_diff(k,w)) then ; 95% confidence interval if (wind_error_diff_avg(k,w).lt.0) then ; then exp is better mres@gsMarkerColor = colors_plot(k+1) else ; ctl is better mres@gsMarkerColor = colors_plot(0) end if gsn_polymarker(wks,plot0,x_axis(w),2.0,mres) end if end if end do end do ;;;;; resbot = True resbot@vpXF = 0.15 ; The left side of the box location resbot@vpYF = 0.20 ; The top side of the plot box location resbot@vpWidthF = 0.70 ; The Width of the plot box resbot@vpHeightF = 0.05 ; The height of the plot box resbot@trYMaxF = 3 ; random value does not need to be changed resbot@trXMinF = 0 resbot@trXMaxF = 28 resbot@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour resbot@tmXBValues = ispan(0,28,2) resbot@tmXBLabels = ispan(0,168,12) resbot@tiYAxisString = "# fcsts" resbot@tiXAxisString = "Forecast lead time (hr)" resbot@tiXAxisFontHeightF = 0.015 resbot@tmYROn = False resbot@tmYLBorderOn = False resbot@tmYLMinorPerMajor = 1 resbot@tmYLOn = False resbot@tmYRBorderOn = False resbot@tmYROn = False resbot@tmYROn = False resbot@tmYRMinorOn = False resbot@tiYAxisFontHeightF = 0.015 resbot@tmXBLabelFontHeightF = 0.015 resbot@tmYLLabelFontHeightF = 0.015 resbot@xyLineThicknessF = 0.05 resbot@tiMainFontHeightF = 0.021 resbot@tmXTOn = False ; turn off the top tickmarks resbot@tmYLOn = False resbot@tmXBMinorOn = False ; No minor tick marks. resbot@gsnDraw = False ; Don't draw individual plot. resbot@gsnFrame = False ; Don't advance frame. ;resbot@gsnYRefLine = 0.0 ; create a reference line ;resbot@gsnAboveYRefLineColor = "green" ; above ref line fill green resbot@gsnXYBarChart = False ; turn on bar chart plot_blank = gsn_csm_xy(wks,x_axis,0.0*fcst_hr_array_ctl(location_max,0:28),resbot) ; plot false values to create plot, then overlay text values on plot ; attache the number of paired cases over the plot y_value_text = (/1.5,1.5,1.5/) do k=0,dimsizes(EXP_name)-1 do w=0,28,2 ; loop over forecast hours, skip first hour for simplicity of viewing image qres = True qres@txFontHeightF = 0.017 qres@txJust = "CenterCenter" qres@txFontColor = colors_plot(k+1) ; color coded by experiment text = gsn_add_text(wks,plot_blank,str_squeeze(tostring(tointeger(fcst_hr_array_exp(k,w)))),x_axis(w),y_value_text(k),qres) end do end do ;;;;;;;;;;;;; draw(plot0) draw(plot_blank) frame(wks) delete(res) delete(bres) delete(resbot) ; intensity wks = gsn_open_wks (output_type,"intensity_error_fcst_hour-all-storms_" + region_file) res = True ; plot mods desired res@gsnXYBarChart = False res@tiMainString = "Intensity Error (" + count_storms + " TCs)" + " " + region_name res@trYMaxF = max_int + 20.0 res@trYMinF = 0 res@gsnFrame = False ; don't advance frame yet res@gsnDraw = False res@gsnFrame = False res@gsnScale = True res@trXMinF = 0 res@trXMaxF = 28 res@vpXF = 0.15 res@vpYF = 0.90 res@vpWidthF = 0.7 res@vpHeightF = 0.6 x_axis = ispan(0,28,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour res@tmXBValues = ispan(0,28,2) res@tmXBLabels = ispan(0,168,12) res@tiYAxisString = "Intensity Error (hPa)" res@tiXAxisString = "Forecast lead time (hr)" res@tiXAxisFontHeightF = 0.015 res@tmYROn = False res@tmYROn = False res@tmYRMinorOn = False res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmXMajorGrid = True ; implement x grid res@tmXMajorGridThicknessF = 1.0 ; 2.0 is default res@tmXMajorGridLineDashPattern = 2 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 8.0 res@tiMainFontHeightF = 0.021 res@xyLineColor = colors_plot(0) plot0 = gsn_csm_xy(wks,x_axis,intensity_error_ctl_fh(location_max,0:28),res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 8.0 plot1 = gsn_csm_xy(wks,x_axis,intensity_error_exp_fh(0,0:28),res) if (num_exp_cases.ge.2) then res@xyLineColor = colors_plot(2) res@xyLineThicknessF = 8.0 plot2 = gsn_csm_xy(wks,x_axis,intensity_error_exp_fh(1,0:28),res) overlay(plot0,plot2) end if if (num_exp_cases.ge.3) then res@xyLineColor = colors_plot(3) res@xyLineThicknessF = 8.0 plot3 = gsn_csm_xy(wks,x_axis,intensity_error_exp_fh(2,0:28),res) overlay(plot0,plot3) end if if (num_exp_cases.ge.4) then res@xyLineColor = colors_plot(4) res@xyLineThicknessF = 8.0 plot4 = gsn_csm_xy(wks,x_axis,intensity_error_exp_fh(3,0:28),res) overlay(plot0,plot4) end if bres = True bres@tmXBMode = "Explicit" bres@gsnFrame = False bres@gsnDraw = False bres@gsnXYBarChart = True bres@gsnXYBarChartColors = "grey78" bres@tmYROn = True bres@tmYRMinorOn = True bres@tmYLOn = False bres@vpXF = 0.15 bres@vpYF = 0.90 bres@vpWidthF = 0.7 bres@vpHeightF = 0.6 bres@trYMinF = 0 bres@tiMainString = "" bres@trYMaxF = max_fcst + 4.0 ; number of samples in each forecast hour bres@tmYRLabelsOn = True bres@tiYAxisFontHeightF = 0.015 bres@tmYRLabelFontHeightF = 0.015 bres@gsnXYBarChartOutlineThicknessF = 1.5 bres@tmYROn = True bres@gsnXYBarChartBarWidth = 0.5 ; change bar widths 0.75 originally bres@tmXBOn = False bres@trXMinF = 0 ; adds space on either end bres@trXMaxF = 28 ; of the 1st and last bars bres@tiYAxisString = "# fcsts" bres@tiYAxisSide = "Right" bres@tiYAxisAngleF = 270 bres@tmXTOn = False plot5 = gsn_csm_xy(wks,x_axis,fcst_hr_array_ctl(location_max,0:28),bres) bres@gsnXYBarChartColors = "grey78" ;plot6 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp(0,0:28),bres) ;draw(plot5) ;draw(plot6) if (num_exp_cases.ge.2) then bres@gsnXYBarChartColors = "grey78" ;plot7 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp(1,0:28),bres) ;draw(plot7) end if if (num_exp_cases.ge.3) then bres@gsnXYBarChartColors = "grey78" ;plot8 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp(2,0:28),bres) ;draw(plot8) end if if (num_exp_cases.ge.4) then bres@gsnXYBarChartColors = "grey78" ;plot9 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp(3,0:28),bres) ;draw(plot9) end if overlay(plot0,plot1) ; Attach a legend marker_value = (/3,4,6,9/) mres = True txres = True txres@txFontHeightF = 0.016 do k=0,dimsizes(labels2)-1 txres@txFontColor = colors_plot(k) gsn_text_ndc(wks,labels2(k),xtxt3(k),ytxt3(k),txres) if (k.ge.1) then mres@gsMarkerIndex = marker_value(k-1) mres@gsMarkerSizeF = 16.0 mres@gsMarkerThicknessF = 4.0 mres@gsMarkerColor = colors_plot(k) gsn_polymarker_ndc(wks,xtxt3(k)-0.07,ytxt3(k),mres) end if end do ; attach significance do k=0,dimsizes(EXP_name)-1 do w=0,28 ; loop over forecast hours mres@gsMarkerIndex = marker_value(k) mres@gsMarkerSizeF = 16.0 mres@gsMarkerThicknessF = 4.0 ;mres@gsMarkerColor = colors_plot(k+1) if (.not.ismissing(intensity_error_diff_avg(k,w)) .and. interval_intensity_diff(k,w).gt.0) then if (abs(intensity_error_diff_avg(k,w)).gt.interval_intensity_diff(k,w)) then ; 95% confidence interval if (intensity_error_diff_avg(k,w).lt.0) then ; then exp is better mres@gsMarkerColor = colors_plot(k+1) else ; ctl is better mres@gsMarkerColor = colors_plot(0) end if gsn_polymarker(wks,plot0,x_axis(w),2.0,mres) end if end if end do end do ;;;;; resbot = True resbot@vpXF = 0.15 ; The left side of the box location resbot@vpYF = 0.20 ; The top side of the plot box location resbot@vpWidthF = 0.70 ; The Width of the plot box resbot@vpHeightF = 0.05 ; The height of the plot box resbot@trYMaxF = 3 ; random value does not need to be changed resbot@trXMinF = 0 resbot@trXMaxF = 28 resbot@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour resbot@tmXBValues = ispan(0,28,2) resbot@tmXBLabels = ispan(0,168,12) resbot@tiYAxisString = "# fcsts" resbot@tiXAxisString = "Forecast lead time (hr)" resbot@tiXAxisFontHeightF = 0.015 resbot@tmYROn = False resbot@tmYROn = False resbot@tmYRMinorOn = False resbot@tmYLBorderOn = False resbot@tmYLMinorPerMajor = 1 resbot@tmYLOn = False resbot@tmYRBorderOn = False resbot@tmYROn = False resbot@tiYAxisFontHeightF = 0.015 resbot@tmXBLabelFontHeightF = 0.015 resbot@tmYLLabelFontHeightF = 0.015 resbot@xyLineThicknessF = 0.05 resbot@tiMainFontHeightF = 0.021 resbot@tmXTOn = False ; turn off the top tickmarks resbot@tmYLOn = False resbot@tmXBMinorOn = False ; No minor tick marks. resbot@gsnDraw = False ; Don't draw individual plot. resbot@gsnFrame = False ; Don't advance frame. ;resbot@gsnYRefLine = 0.0 ; create a reference line ;resbot@gsnAboveYRefLineColor = "green" ; above ref line fill green resbot@gsnXYBarChart = False ; turn on bar chart plot_blank = gsn_csm_xy(wks,x_axis,0.0*fcst_hr_array_ctl(location_max,0:28),resbot) ; plot false values to create plot, then overlay text values on plot ; attache the number of paired cases over the plot y_value_text = (/1.5,1.5,1.5/) do k=0,dimsizes(EXP_name)-1 do w=0,28,2 ; loop over forecast hours, skip first hour for simplicity of viewing image qres = True qres@txFontHeightF = 0.017 qres@txJust = "CenterCenter" qres@txFontColor = colors_plot(k+1) ; color coded by experiment text = gsn_add_text(wks,plot_blank,str_squeeze(tostring(tointeger(fcst_hr_array_exp(k,w)))),x_axis(w),y_value_text(k),qres) end do end do ;;;;;;;;;;;;; draw(plot0) draw(plot_blank) frame(wks) delete(res) delete(bres) delete(resbot) ; do the same 3 plots, but include instead the 1-sigma standard deviation for uncertainty wks = gsn_open_wks (output_type,"track_error_fcst_hour-all-storms_stddev_" + region_file) res = True ; plot mods desired res@gsnXYBarChart = False res@tiMainString = "Track Error (" + count_storms + " TCs)" + " " + region_name res@trYMaxF = max_track + 200.0 res@trYMinF = 0 res@gsnFrame = False ; don't advance frame yet res@gsnDraw = False res@gsnFrame = False res@vpXF = 0.15 res@vpYF = 0.90 res@vpWidthF = 0.7 res@vpHeightF = 0.6 res@gsnScale = True res@trXMinF = 0 res@trXMaxF = 28 x_axis = ispan(0,28,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour res@tmXBValues = ispan(0,28,2) res@tmXBLabels = ispan(0,168,12) res@tiYAxisString = "Track Error (km)" res@tiXAxisString = "Forecast lead time (hr)" res@tiXAxisFontHeightF = 0.015 res@tmYROn = False res@tmYROn = False res@tmYRMinorOn = False res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmXMajorGrid = True ; implement x grid res@tmXMajorGridThicknessF = 1.0 ; 2.0 is default res@tmXMajorGridLineDashPattern = 2 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 8.0 res@tiMainFontHeightF = 0.021 res@xyLineColor = colors_plot(0) plot0 = gsn_csm_xy(wks,x_axis,track_error_ctl_fh(location_max,0:28),res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 8.0 plot1 = gsn_csm_xy(wks,x_axis,track_error_exp_fh(0,0:28),res) overlay(plot0,plot1) if (num_exp_cases.ge.2) then res@xyLineColor = colors_plot(2) res@xyLineThicknessF = 8.0 plot2 = gsn_csm_xy(wks,x_axis,track_error_exp_fh(1,0:28),res) overlay(plot0,plot2) end if if (num_exp_cases.ge.3) then res@xyLineColor = colors_plot(3) res@xyLineThicknessF = 8.0 plot3 = gsn_csm_xy(wks,x_axis,track_error_exp_fh(2,0:28),res) overlay(plot0,plot3) end if if (num_exp_cases.ge.4) then res@xyLineColor = colors_plot(4) res@xyLineThicknessF = 8.0 plot4 = gsn_csm_xy(wks,x_axis,track_error_exp_fh(3,0:28),res) overlay(plot0,plot4) end if ; add standard deviation spread res@gsnXYFillColors = "gray78" res@xyLineColor = -1 ind_valid = ind(.not.ismissing(track_error_ctl_fh_std(location_max,:))) max_std = (track_error_ctl_fh_std(location_max,ind_valid)) / sqrt(fcst_hr_array_ctl(location_max,ind_valid)) spread = (/track_error_ctl_fh(location_max,ind_valid)-max_std,track_error_ctl_fh(location_max,ind_valid)+max_std/) res@gsnXYFillOpacities = 0.45 ctl_spread = gsn_csm_xy(wks,x_axis(ind_valid),spread,res) res@gsnXYFillColors = colors_plot(1) res@xyLineColor = -1 delete(ind_valid) delete(max_std) delete(spread) ind_valid = ind(.not.ismissing(track_error_exp_fh_std(0,:))) max_std = (track_error_exp_fh_std(0,ind_valid)) / sqrt(fcst_hr_array_exp(0,ind_valid)) spread = (/track_error_exp_fh(0,ind_valid)-max_std,track_error_exp_fh(0,ind_valid)+max_std/) res@gsnXYFillOpacities = 0.2 exp_spread1 = gsn_csm_xy(wks,x_axis(ind_valid),spread,res) overlay(plot0,ctl_spread) overlay(plot0,exp_spread1) if (num_exp_cases.ge.2) then res@gsnXYFillColors = colors_plot(2) res@xyLineColor = -1 delete(ind_valid) delete(max_std) delete(spread) ind_valid = ind(.not.ismissing(track_error_exp_fh_std(1,:))) max_std = (track_error_exp_fh_std(1,ind_valid)) / sqrt(fcst_hr_array_exp(1,ind_valid)) spread = (/track_error_exp_fh(1,ind_valid)-max_std,track_error_exp_fh(1,ind_valid)+max_std/) res@gsnXYFillOpacities = 0.2 exp_spread2 = gsn_csm_xy(wks,x_axis(ind_valid),spread,res) overlay(plot0,exp_spread2) end if if (num_exp_cases.ge.3) then res@gsnXYFillColors = colors_plot(3) res@xyLineColor = -1 delete(ind_valid) delete(max_std) delete(spread) ind_valid = ind(.not.ismissing(track_error_exp_fh_std(2,:))) max_std = (track_error_exp_fh_std(2,ind_valid)) / sqrt(fcst_hr_array_exp(2,ind_valid)) spread = (/track_error_exp_fh(2,ind_valid)-max_std,track_error_exp_fh(2,ind_valid)+max_std/) res@gsnXYFillOpacities = 0.2 exp_spread3 = gsn_csm_xy(wks,x_axis(ind_valid),spread,res) overlay(plot0,exp_spread3) end if if (num_exp_cases.ge.4) then res@gsnXYFillColors = colors_plot(4) res@xyLineColor = -1 delete(ind_valid) delete(max_std) delete(spread) ind_valid = ind(.not.ismissing(track_error_exp_fh_std(3,:))) max_std = (track_error_exp_fh_std(3,ind_valid)) / sqrt(fcst_hr_array_exp(3,ind_valid)) spread = (/track_error_exp_fh(3,ind_valid)-max_std,track_error_exp_fh(3,ind_valid)+max_std/) res@gsnXYFillOpacities = 0.2 exp_spread4 = gsn_csm_xy(wks,x_axis(ind_valid),spread,res) overlay(plot0,exp_spread4) end if bres = True bres@tmXBMode = "Explicit" bres@gsnFrame = False bres@gsnDraw = False bres@gsnXYBarChart = True bres@gsnXYBarChartColors = colors_plot(0) bres@tmYROn = True bres@tmYRMinorOn = True bres@tmYLOn = False bres@vpXF = 0.15 bres@vpYF = 0.90 bres@vpWidthF = 0.7 bres@vpHeightF = 0.6 bres@trYMinF = 0 bres@tiMainString = "" bres@trYMaxF = max_fcst + 4.0 ; number of samples in each forecast hour bres@tmYRLabelsOn = True bres@tiYAxisFontHeightF = 0.015 bres@tmYRLabelFontHeightF = 0.015 bres@gsnXYBarChartOutlineThicknessF = 1.5 bres@tmYROn = True bres@gsnXYBarChartBarWidth = 0.5 ; change bar widths 0.75 originally bres@tmXBOn = False bres@trXMinF = 0 ; adds space on either end bres@trXMaxF = 28 ; of the 1st and last bars bres@tiYAxisString = "# fcsts" bres@tiYAxisSide = "Right" bres@tiYAxisAngleF = 270 bres@tmXTOn = False bres@gsnXYBarChartFillOpacityF = 0.5 ;plot5 = gsn_csm_xy(wks,x_axis,fcst_hr_array_ctl(location_max,0:28),bres) bres@gsnXYBarChartColors = colors_plot(1) ;plot6 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp(0,0:28),bres) ;draw(plot5) ;draw(plot6) if (num_exp_cases.ge.2) then bres@gsnXYBarChartColors = colors_plot(2) ;plot7 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp(1,0:28),bres) ;draw(plot7) end if if (num_exp_cases.ge.3) then bres@gsnXYBarChartColors = colors_plot(3) ;plot8 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp(2,0:28),bres) ;draw(plot8) end if if (num_exp_cases.ge.4) then bres@gsnXYBarChartColors = colors_plot(4) ;plot9 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp(3,0:28),bres) ;draw(plot9) end if ; Attach a legend marker_value = (/3,4,6,9/) mres = True txres = True txres@txFontHeightF = 0.016 do k=0,dimsizes(labels2)-1 txres@txFontColor = colors_plot(k) gsn_text_ndc(wks,labels2(k),xtxt3(k),ytxt3(k),txres) if (k.ge.1) then mres@gsMarkerIndex = marker_value(k-1) mres@gsMarkerSizeF = 16.0 mres@gsMarkerThicknessF = 4.0 mres@gsMarkerColor = colors_plot(k) gsn_polymarker_ndc(wks,xtxt3(k)-0.07,ytxt3(k),mres) end if end do ; attach significance do k=0,dimsizes(EXP_name)-1 do w=0,28 ; loop over forecast hours mres@gsMarkerIndex = marker_value(k) mres@gsMarkerSizeF = 16.0 mres@gsMarkerThicknessF = 4.0 ;mres@gsMarkerColor = colors_plot(k+1) if (.not.ismissing(track_error_diff_avg(k,w)) .and. interval_track_diff(k,w).gt.0) then if (abs(track_error_diff_avg(k,w)).gt.interval_track_diff(k,w)) then ; 95% confidence interval if (track_error_diff_avg(k,w).lt.0) then ; then exp is better mres@gsMarkerColor = colors_plot(k+1) else ; ctl is better mres@gsMarkerColor = colors_plot(0) end if gsn_polymarker(wks,plot0,x_axis(w),20.0,mres) end if end if end do end do ;;;;; resbot = True resbot@vpXF = 0.15 ; The left side of the box location resbot@vpYF = 0.20 ; The top side of the plot box location resbot@vpWidthF = 0.70 ; The Width of the plot box resbot@vpHeightF = 0.05 ; The height of the plot box resbot@trYMaxF = 3 ; random value does not need to be changed resbot@trXMinF = 0 resbot@trXMaxF = 28 resbot@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour resbot@tmXBValues = ispan(0,28,2) resbot@tmXBLabels = ispan(0,168,12) resbot@tiYAxisString = "# fcsts" resbot@tiXAxisString = "Forecast lead time (hr)" resbot@tiXAxisFontHeightF = 0.015 resbot@tmYROn = False resbot@tmYLBorderOn = False resbot@tmYLMinorPerMajor = 1 resbot@tmYLOn = False resbot@tmYRBorderOn = False resbot@tmYROn = False resbot@tmYROn = False resbot@tmYRMinorOn = False resbot@tiYAxisFontHeightF = 0.015 resbot@tmXBLabelFontHeightF = 0.015 resbot@tmYLLabelFontHeightF = 0.015 resbot@xyLineThicknessF = 0.05 resbot@tiMainFontHeightF = 0.021 resbot@tmXTOn = False ; turn off the top tickmarks resbot@tmYLOn = False resbot@tmXBMinorOn = False ; No minor tick marks. resbot@gsnDraw = False ; Don't draw individual plot. resbot@gsnFrame = False ; Don't advance frame. ;resbot@gsnYRefLine = 0.0 ; create a reference line ;resbot@gsnAboveYRefLineColor = "green" ; above ref line fill green resbot@gsnXYBarChart = False ; turn on bar chart plot_blank = gsn_csm_xy(wks,x_axis,0.0*fcst_hr_array_ctl(location_max,0:28),resbot) ; plot false values to create plot, then overlay text values on plot ; attache the number of paired cases over the plot y_value_text = (/1.5,1.5,1.5/) do k=0,dimsizes(EXP_name)-1 do w=0,28,2 ; loop over forecast hours, skip first hour for simplicity of viewing image qres = True qres@txFontHeightF = 0.017 qres@txJust = "CenterCenter" qres@txFontColor = colors_plot(k+1) ; color coded by experiment text = gsn_add_text(wks,plot_blank,str_squeeze(tostring(tointeger(fcst_hr_array_exp(k,w)))),x_axis(w),y_value_text(k),qres) end do end do ;;;;;;;;;;;;; draw(plot0) draw(plot_blank) frame(wks) delete(res) delete(bres) delete(resbot) ; wind wks = gsn_open_wks (output_type,"wind_error_fcst_hour-all-storms_stddev_" + region_file) res = True ; plot mods desired res@gsnXYBarChart = False res@tiMainString = "Wind Error (" + count_storms + " TCs)" + " " + region_name res@trYMaxF = max_wind + 20.0 res@trYMinF = 0 res@gsnFrame = False ; don't advance frame yet res@gsnDraw = False res@gsnFrame = False res@vpXF = 0.15 res@vpYF = 0.90 res@vpWidthF = 0.7 res@vpHeightF = 0.6 res@gsnScale = True res@trXMinF = 0 res@trXMaxF = 28 x_axis = ispan(0,28,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour res@tmXBValues = ispan(0,28,2) res@tmXBLabels = ispan(0,168,12) res@tiYAxisString = "Wind Error (kts)" res@tiXAxisString = "Forecast lead time (hr)" res@tiXAxisFontHeightF = 0.015 res@tmYROn = False res@tmYROn = False res@tmYRMinorOn = False res@tiYAxisFontHeightF = 0.015 res@tmXMajorGrid = True ; implement x grid res@tmXMajorGridThicknessF = 1.0 ; 2.0 is default res@tmXMajorGridLineDashPattern = 2 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 8.0 res@tiMainFontHeightF = 0.021 res@xyLineColor = colors_plot(0) plot0 = gsn_csm_xy(wks,x_axis,wind_error_ctl_fh(location_max,0:28),res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 8.0 plot1 = gsn_csm_xy(wks,x_axis,wind_error_exp_fh(0,0:28),res) overlay(plot0,plot1) if (num_exp_cases.ge.2) then res@xyLineColor = colors_plot(2) res@xyLineThicknessF = 8.0 plot2 = gsn_csm_xy(wks,x_axis,wind_error_exp_fh(1,0:28),res) overlay(plot0,plot2) end if if (num_exp_cases.ge.3) then res@xyLineColor = colors_plot(3) res@xyLineThicknessF = 8.0 plot3 = gsn_csm_xy(wks,x_axis,wind_error_exp_fh(2,0:28),res) overlay(plot0,plot3) end if if (num_exp_cases.ge.4) then res@xyLineColor = colors_plot(4) res@xyLineThicknessF = 8.0 plot4 = gsn_csm_xy(wks,x_axis,wind_error_exp_fh(3,0:28),res) overlay(plot0,plot4) end if ; add standard deviation spread res@gsnXYFillColors = "gray78" res@xyLineColor = -1 delete(ind_valid) delete(max_std) delete(spread) ind_valid = ind(.not.ismissing(wind_error_ctl_fh_std(location_max,:))) max_std = (wind_error_ctl_fh_std(location_max,ind_valid)) / sqrt(fcst_hr_array_ctl(location_max,ind_valid)) spread = (/wind_error_ctl_fh(location_max,ind_valid)-max_std,wind_error_ctl_fh(location_max,ind_valid)+max_std/) res@gsnXYFillOpacities = 0.45 ctl_spread = gsn_csm_xy(wks,x_axis(ind_valid),spread,res) res@gsnXYFillColors = colors_plot(1) res@xyLineColor = -1 delete(ind_valid) delete(max_std) delete(spread) ind_valid = ind(.not.ismissing(wind_error_exp_fh_std(0,:))) max_std = (wind_error_exp_fh_std(0,ind_valid)) / sqrt(fcst_hr_array_exp(0,ind_valid)) spread = (/wind_error_exp_fh(0,ind_valid)-max_std,wind_error_exp_fh(0,ind_valid)+max_std/) res@gsnXYFillOpacities = 0.2 exp_spread1 = gsn_csm_xy(wks,x_axis(ind_valid),spread,res) overlay(plot0,ctl_spread) overlay(plot0,exp_spread1) if (num_exp_cases.ge.2) then res@gsnXYFillColors = colors_plot(2) res@xyLineColor = -1 delete(ind_valid) delete(max_std) delete(spread) ind_valid = ind(.not.ismissing(wind_error_exp_fh_std(1,:))) max_std = (wind_error_exp_fh_std(1,ind_valid)) / sqrt(fcst_hr_array_exp(1,ind_valid)) spread = (/wind_error_exp_fh(1,ind_valid)-max_std,wind_error_exp_fh(1,ind_valid)+max_std/) res@gsnXYFillOpacities = 0.2 exp_spread2 = gsn_csm_xy(wks,x_axis(ind_valid),spread,res) overlay(plot0,exp_spread2) end if if (num_exp_cases.ge.3) then res@gsnXYFillColors = colors_plot(3) res@xyLineColor = -1 delete(ind_valid) delete(max_std) delete(spread) ind_valid = ind(.not.ismissing(wind_error_exp_fh_std(2,:))) max_std = (wind_error_exp_fh_std(2,ind_valid)) / sqrt(fcst_hr_array_exp(2,ind_valid)) spread = (/wind_error_exp_fh(2,ind_valid)-max_std,wind_error_exp_fh(2,ind_valid)+max_std/) res@gsnXYFillOpacities = 0.2 exp_spread3 = gsn_csm_xy(wks,x_axis(ind_valid),spread,res) overlay(plot0,exp_spread3) end if if (num_exp_cases.ge.4) then res@gsnXYFillColors = colors_plot(4) res@xyLineColor = -1 delete(ind_valid) delete(max_std) delete(spread) ind_valid = ind(.not.ismissing(wind_error_exp_fh_std(3,:))) max_std = (wind_error_exp_fh_std(3,ind_valid)) / sqrt(fcst_hr_array_exp(3,ind_valid)) spread = (/wind_error_exp_fh(3,ind_valid)-max_std,wind_error_exp_fh(3,ind_valid)+max_std/) res@gsnXYFillOpacities = 0.2 exp_spread4 = gsn_csm_xy(wks,x_axis(ind_valid),spread,res) overlay(plot0,exp_spread4) end if bres = True bres@tmXBMode = "Explicit" bres@gsnFrame = False bres@gsnDraw = False bres@gsnXYBarChart = True bres@gsnXYBarChartColors = "grey78" bres@tmYROn = True bres@tmYRMinorOn = True bres@tmYLOn = False bres@vpXF = 0.15 bres@vpYF = 0.90 bres@vpWidthF = 0.7 bres@vpHeightF = 0.6 bres@trYMinF = 0 bres@tiMainString = "" bres@trYMaxF = max_fcst + 4.0 ; number of samples in each forecast hour bres@tmYRLabelsOn = True bres@tiYAxisFontHeightF = 0.015 bres@tmYRLabelFontHeightF = 0.015 bres@gsnXYBarChartOutlineThicknessF = 1.5 bres@tmYROn = True bres@gsnXYBarChartBarWidth = 0.5 ; change bar widths 0.75 originally bres@tmXBOn = False bres@trXMinF = 0 ; adds space on either end bres@trXMaxF = 28 ; of the 1st and last bars bres@tiYAxisString = "# fcsts" bres@tiYAxisSide = "Right" bres@tiYAxisAngleF = 270 bres@tmXTOn = False plot5 = gsn_csm_xy(wks,x_axis,fcst_hr_array_ctl(location_max,0:28),bres) bres@gsnXYBarChartColors = "grey78" ;plot6 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp(0,0:28),bres) ;draw(plot5) ;draw(plot6) if (num_exp_cases.ge.2) then bres@gsnXYBarChartColors = "grey78" ;plot7 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp(1,0:28),bres) ;draw(plot7) end if if (num_exp_cases.ge.3) then bres@gsnXYBarChartColors = "grey78" ;plot8 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp(2,0:28),bres) ;draw(plot8) end if if (num_exp_cases.ge.4) then bres@gsnXYBarChartColors = "grey78" ;plot9 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp(3,0:28),bres) ;draw(plot9) end if ; Attach a legend marker_value = (/3,4,6,9/) mres = True txres = True txres@txFontHeightF = 0.016 do k=0,dimsizes(labels2)-1 txres@txFontColor = colors_plot(k) gsn_text_ndc(wks,labels2(k),xtxt3(k),ytxt3(k),txres) if (k.ge.1) then mres@gsMarkerIndex = marker_value(k-1) mres@gsMarkerSizeF = 16.0 mres@gsMarkerThicknessF = 4.0 mres@gsMarkerColor = colors_plot(k) gsn_polymarker_ndc(wks,xtxt3(k)-0.07,ytxt3(k),mres) end if end do ; attach significance do k=0,dimsizes(EXP_name)-1 do w=0,28 ; loop over forecast hours mres@gsMarkerIndex = marker_value(k) mres@gsMarkerSizeF = 16.0 mres@gsMarkerThicknessF = 4.0 ;mres@gsMarkerColor = colors_plot(k+1) if (.not.ismissing(wind_error_diff_avg(k,w)) .and. interval_wind_diff(k,w).gt.0) then if (abs(wind_error_diff_avg(k,w)).gt.interval_wind_diff(k,w)) then ; 95% confidence interval if (wind_error_diff_avg(k,w).lt.0) then ; then exp is better mres@gsMarkerColor = colors_plot(k+1) else ; ctl is better mres@gsMarkerColor = colors_plot(0) end if gsn_polymarker(wks,plot0,x_axis(w),2.0,mres) end if end if end do end do ;;;;; resbot = True resbot@vpXF = 0.15 ; The left side of the box location resbot@vpYF = 0.20 ; The top side of the plot box location resbot@vpWidthF = 0.70 ; The Width of the plot box resbot@vpHeightF = 0.05 ; The height of the plot box resbot@trYMaxF = 3 ; random value does not need to be changed resbot@trXMinF = 0 resbot@trXMaxF = 28 resbot@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour resbot@tmXBValues = ispan(0,28,2) resbot@tmXBLabels = ispan(0,168,12) resbot@tiYAxisString = "# fcsts" resbot@tiXAxisString = "Forecast lead time (hr)" resbot@tiXAxisFontHeightF = 0.015 resbot@tmYLBorderOn = False resbot@tmYLMinorPerMajor = 1 resbot@tmYLOn = False resbot@tmYRBorderOn = False resbot@tmYROn = False resbot@tmYROn = False resbot@tmYROn = False resbot@tmYRMinorOn = False resbot@tiYAxisFontHeightF = 0.015 resbot@tmXBLabelFontHeightF = 0.015 resbot@tmYLLabelFontHeightF = 0.015 resbot@xyLineThicknessF = 0.05 resbot@tiMainFontHeightF = 0.021 resbot@tmXTOn = False ; turn off the top tickmarks resbot@tmYLOn = False resbot@tmXBMinorOn = False ; No minor tick marks. resbot@gsnDraw = False ; Don't draw individual plot. resbot@gsnFrame = False ; Don't advance frame. ;resbot@gsnYRefLine = 0.0 ; create a reference line ;resbot@gsnAboveYRefLineColor = "green" ; above ref line fill green resbot@gsnXYBarChart = False ; turn on bar chart plot_blank = gsn_csm_xy(wks,x_axis,0.0*fcst_hr_array_ctl(location_max,0:28),resbot) ; plot false values to create plot, then overlay text values on plot ; attache the number of paired cases over the plot y_value_text = (/1.5,1.5,1.5/) do k=0,dimsizes(EXP_name)-1 do w=0,28,2 ; loop over forecast hours, skip first hour for simplicity of viewing image qres = True qres@txFontHeightF = 0.017 qres@txJust = "CenterCenter" qres@txFontColor = colors_plot(k+1) ; color coded by experiment text = gsn_add_text(wks,plot_blank,str_squeeze(tostring(tointeger(fcst_hr_array_exp(k,w)))),x_axis(w),y_value_text(k),qres) end do end do ;;;;;;;;;;;;; draw(plot0) draw(plot_blank) frame(wks) delete(res) delete(bres) delete(resbot) ; intensity wks = gsn_open_wks (output_type,"intensity_error_fcst_hour-all-storms_stddev_" + region_file) res = True ; plot mods desired res@gsnXYBarChart = False res@tiMainString = "Intensity Error (" + count_storms + " TCs)" + " " + region_name res@trYMaxF = max_int + 20.0 res@trYMinF = 0 res@gsnFrame = False ; don't advance frame yet res@gsnDraw = False res@gsnFrame = False res@vpXF = 0.15 res@vpYF = 0.90 res@vpWidthF = 0.7 res@vpHeightF = 0.6 res@gsnScale = True res@trXMinF = 0 res@trXMaxF = 28 x_axis = ispan(0,28,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour res@tmXBValues = ispan(0,28,2) res@tmXBLabels = ispan(0,168,12) res@tiYAxisString = "Intensity Error (hPa)" res@tiXAxisString = "Forecast lead time (hr)" res@tiXAxisFontHeightF = 0.015 res@tmYROn = False res@tmYROn = False res@tmXMajorGrid = True ; implement x grid res@tmXMajorGridThicknessF = 1.0 ; 2.0 is default res@tmXMajorGridLineDashPattern = 2 res@tmYRMinorOn = False res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 8.0 res@tiMainFontHeightF = 0.021 res@xyLineColor = colors_plot(0) plot0 = gsn_csm_xy(wks,x_axis,intensity_error_ctl_fh(location_max,0:28),res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 8.0 plot1 = gsn_csm_xy(wks,x_axis,intensity_error_exp_fh(0,0:28),res) overlay(plot0,plot1) if (num_exp_cases.ge.2) then res@xyLineColor = colors_plot(2) res@xyLineThicknessF = 8.0 plot2 = gsn_csm_xy(wks,x_axis,intensity_error_exp_fh(1,0:28),res) overlay(plot0,plot2) end if if (num_exp_cases.ge.3) then res@xyLineColor = colors_plot(3) res@xyLineThicknessF = 8.0 plot3 = gsn_csm_xy(wks,x_axis,intensity_error_exp_fh(2,0:28),res) overlay(plot0,plot3) end if if (num_exp_cases.ge.4) then res@xyLineColor = colors_plot(4) res@xyLineThicknessF = 8.0 plot4 = gsn_csm_xy(wks,x_axis,intensity_error_exp_fh(3,0:28),res) overlay(plot0,plot4) end if ; add standard deviation spread res@gsnXYFillColors = "gray78" res@xyLineColor = -1 delete(ind_valid) delete(max_std) delete(spread) ind_valid = ind(.not.ismissing(intensity_error_ctl_fh_std(location_max,:))) max_std = (intensity_error_ctl_fh_std(location_max,ind_valid)) / sqrt(fcst_hr_array_ctl(location_max,ind_valid)) spread = (/intensity_error_ctl_fh(location_max,ind_valid)-max_std,intensity_error_ctl_fh(location_max,ind_valid)+max_std/) res@gsnXYFillOpacities = 0.45 ctl_spread = gsn_csm_xy(wks,x_axis(ind_valid),spread,res) res@gsnXYFillColors = colors_plot(1) res@xyLineColor = -1 delete(ind_valid) delete(max_std) delete(spread) ind_valid = ind(.not.ismissing(intensity_error_exp_fh_std(0,:))) max_std = (intensity_error_exp_fh_std(0,ind_valid)) / sqrt(fcst_hr_array_exp(0,ind_valid)) spread = (/intensity_error_exp_fh(0,ind_valid)-max_std,intensity_error_exp_fh(0,ind_valid)+max_std/) res@gsnXYFillOpacities = 0.2 exp_spread1 = gsn_csm_xy(wks,x_axis(ind_valid),spread,res) overlay(plot0,ctl_spread) overlay(plot0,exp_spread1) if (num_exp_cases.ge.2) then res@gsnXYFillColors = colors_plot(2) res@xyLineColor = -1 delete(ind_valid) delete(max_std) delete(spread) ind_valid = ind(.not.ismissing(intensity_error_exp_fh_std(1,:))) max_std = (intensity_error_exp_fh_std(1,ind_valid)) / sqrt(fcst_hr_array_exp(1,ind_valid)) spread = (/intensity_error_exp_fh(1,ind_valid)-max_std,intensity_error_exp_fh(1,ind_valid)+max_std/) res@gsnXYFillOpacities = 0.2 exp_spread2 = gsn_csm_xy(wks,x_axis(ind_valid),spread,res) overlay(plot0,exp_spread2) end if if (num_exp_cases.ge.3) then res@gsnXYFillColors = colors_plot(3) res@xyLineColor = -1 delete(ind_valid) delete(max_std) delete(spread) ind_valid = ind(.not.ismissing(intensity_error_exp_fh_std(2,:))) max_std = (intensity_error_exp_fh_std(2,ind_valid)) / sqrt(fcst_hr_array_exp(2,ind_valid)) spread = (/intensity_error_exp_fh(2,ind_valid)-max_std,intensity_error_exp_fh(2,ind_valid)+max_std/) res@gsnXYFillOpacities = 0.2 exp_spread3 = gsn_csm_xy(wks,x_axis(ind_valid),spread,res) overlay(plot0,exp_spread3) end if if (num_exp_cases.ge.4) then res@gsnXYFillColors = colors_plot(4) res@xyLineColor = -1 delete(ind_valid) delete(max_std) delete(spread) ind_valid = ind(.not.ismissing(intensity_error_exp_fh_std(3,:))) max_std = (intensity_error_exp_fh_std(3,ind_valid)) / sqrt(fcst_hr_array_exp(3,ind_valid)) spread = (/intensity_error_exp_fh(3,ind_valid)-max_std,intensity_error_exp_fh(3,ind_valid)+max_std/) res@gsnXYFillOpacities = 0.2 exp_spread4 = gsn_csm_xy(wks,x_axis(ind_valid),spread,res) overlay(plot0,exp_spread4) end if bres = True bres@tmXBMode = "Explicit" bres@gsnFrame = False bres@gsnDraw = False bres@gsnXYBarChart = True bres@gsnXYBarChartColors = "grey78" bres@tmYROn = True bres@tmYRMinorOn = True bres@vpXF = 0.15 bres@vpYF = 0.90 bres@vpWidthF = 0.7 bres@vpHeightF = 0.6 bres@tmYLOn = False bres@trYMinF = 0 bres@tiMainString = "" bres@trYMaxF = max_fcst + 4.0 ; number of samples in each forecast hour bres@tmYRLabelsOn = True bres@tiYAxisFontHeightF = 0.015 bres@tmYRLabelFontHeightF = 0.015 bres@gsnXYBarChartOutlineThicknessF = 1.5 bres@tmYROn = True bres@gsnXYBarChartBarWidth = 0.5 ; change bar widths 0.75 originally bres@tmXBOn = False bres@trXMinF = 0 ; adds space on either end bres@trXMaxF = 28 ; of the 1st and last bars bres@tiYAxisString = "# fcsts" bres@tiYAxisSide = "Right" bres@tiYAxisAngleF = 270 bres@tmXTOn = False plot5 = gsn_csm_xy(wks,x_axis,fcst_hr_array_ctl(location_max,0:28),bres) bres@gsnXYBarChartColors = "grey78" ;plot6 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp(0,0:28),bres) ;draw(plot5) ;draw(plot6) if (num_exp_cases.ge.2) then bres@gsnXYBarChartColors = "grey78" ;plot7 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp(1,0:28),bres) ;draw(plot7) end if if (num_exp_cases.ge.3) then bres@gsnXYBarChartColors = "grey78" ;plot8 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp(2,0:28),bres) ;draw(plot8) end if if (num_exp_cases.ge.4) then bres@gsnXYBarChartColors = "grey78" ;plot9 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp(3,0:28),bres) ;draw(plot9) end if ; Attach a legend marker_value = (/3,4,6,9/) mres = True txres = True txres@txFontHeightF = 0.016 do k=0,dimsizes(labels2)-1 txres@txFontColor = colors_plot(k) gsn_text_ndc(wks,labels2(k),xtxt3(k),ytxt3(k),txres) if (k.ge.1) then mres@gsMarkerIndex = marker_value(k-1) mres@gsMarkerSizeF = 16.0 mres@gsMarkerThicknessF = 4.0 mres@gsMarkerColor = colors_plot(k) gsn_polymarker_ndc(wks,xtxt3(k)-0.07,ytxt3(k),mres) end if end do ; attach significance do k=0,dimsizes(EXP_name)-1 do w=0,28 ; loop over forecast hours mres@gsMarkerIndex = marker_value(k) mres@gsMarkerSizeF = 16.0 mres@gsMarkerThicknessF = 4.0 ;mres@gsMarkerColor = colors_plot(k+1) if (.not.ismissing(intensity_error_diff_avg(k,w)) .and. interval_intensity_diff(k,w).gt.0) then if (abs(intensity_error_diff_avg(k,w)).gt.interval_intensity_diff(k,w)) then ; 95% confidence interval if (intensity_error_diff_avg(k,w).lt.0) then ; then exp is better mres@gsMarkerColor = colors_plot(k+1) else ; ctl is better mres@gsMarkerColor = colors_plot(0) end if gsn_polymarker(wks,plot0,x_axis(w),2.0,mres) end if end if end do end do ;;;;; resbot = True resbot@vpXF = 0.15 ; The left side of the box location resbot@vpYF = 0.20 ; The top side of the plot box location resbot@vpWidthF = 0.70 ; The Width of the plot box resbot@vpHeightF = 0.05 ; The height of the plot box resbot@trYMaxF = 3 ; random value does not need to be changed resbot@trXMinF = 0 resbot@trXMaxF = 28 resbot@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour resbot@tmXBValues = ispan(0,28,2) resbot@tmXBLabels = ispan(0,168,12) resbot@tiYAxisString = "# fcsts" resbot@tiXAxisString = "Forecast lead time (hr)" resbot@tiXAxisFontHeightF = 0.015 resbot@tmYLBorderOn = False resbot@tmYLMinorPerMajor = 1 resbot@tmYLOn = False resbot@tmYRBorderOn = False resbot@tmYROn = False resbot@tmYROn = False resbot@tmYROn = False resbot@tmYRMinorOn = False resbot@tiYAxisFontHeightF = 0.015 resbot@tmXBLabelFontHeightF = 0.015 resbot@tmYLLabelFontHeightF = 0.015 resbot@xyLineThicknessF = 0.05 resbot@tiMainFontHeightF = 0.021 resbot@tmXTOn = False ; turn off the top tickmarks resbot@tmYLOn = False resbot@tmXBMinorOn = False ; No minor tick marks. resbot@gsnDraw = False ; Don't draw individual plot. resbot@gsnFrame = False ; Don't advance frame. ;resbot@gsnYRefLine = 0.0 ; create a reference line ;resbot@gsnAboveYRefLineColor = "green" ; above ref line fill green resbot@gsnXYBarChart = False ; turn on bar chart plot_blank = gsn_csm_xy(wks,x_axis,0.0*fcst_hr_array_ctl(location_max,0:28),resbot) ; plot false values to create plot, then overlay text values on plot ; attache the number of paired cases over the plot y_value_text = (/1.5,1.5,1.5/) do k=0,dimsizes(EXP_name)-1 do w=0,28,2 ; loop over forecast hours, skip first hour for simplicity of viewing image qres = True qres@txFontHeightF = 0.017 qres@txJust = "CenterCenter" qres@txFontColor = colors_plot(k+1) ; color coded by experiment text = gsn_add_text(wks,plot_blank,str_squeeze(tostring(tointeger(fcst_hr_array_exp(k,w)))),x_axis(w),y_value_text(k),qres) end do end do ;;;;;;;;;;;;; draw(plot0) draw(plot_blank) frame(wks) delete(res) delete(bres) delete(resbot) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; PLOT AVERAGE DIFFERENCE AND 95% CONFIDENCE INTERVAL FOR TRACK, WIND, AND INTENSITY ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; delete(ind_valid) ; intensity error max_diff = max(intensity_error_diff_avg) ; set bounds min_diff = min(intensity_error_diff_avg) range_diff = max_diff - min_diff wks = gsn_open_wks (output_type,"intensity_confidence_interval_fcst_hour-all-storms_" + region_file) res = True ; plot mods desired res@gsnXYBarChart = False res@tiMainString = "Intensity Error Difference (" + count_storms + " TCs)" + " " + region_name res@trYMaxF = max_diff + 1.2*range_diff res@trYMinF = min_diff - 1.2*range_diff res@gsnFrame = False ; don't advance frame yet res@gsnDraw = False res@gsnFrame = False res@vpXF = 0.15 res@vpYF = 0.90 res@vpWidthF = 0.7 res@vpHeightF = 0.6 res@gsnScale = True res@trXMinF = 0 res@trXMaxF = 28 x_axis = ispan(0,28,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour res@tmXBValues = ispan(0,28,2) res@tmXBLabels = ispan(0,168,12) res@tiYAxisString = "Intensity Error Difference (hPa)" res@tiXAxisString = "Forecast lead time (hr)" res@tiXAxisFontHeightF = 0.015 res@tmYROn = False res@tmYROn = False res@tmXMajorGrid = True ; implement x grid res@tmXMajorGridThicknessF = 1.0 ; 2.0 is default res@tmXMajorGridLineDashPattern = 2 res@tmYRMinorOn = False res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 8.0 res@tiMainFontHeightF = 0.021 res@xyLineColor = colors_plot(1) res@gsnYRefLine = (/0/) ; first plot ind_valid = ind(.not.ismissing(intensity_error_diff_avg(0,:))) plot0 = gsn_csm_xy(wks,x_axis(ind_valid),intensity_error_diff_avg(0,ind_valid),res) if (num_exp_cases.ge.2) then delete(ind_valid) ind_valid = ind(.not.ismissing(intensity_error_diff_avg(1,:))) res@xyLineColor = colors_plot(2) res@xyLineThicknessF = 8.0 plot2 = gsn_csm_xy(wks,x_axis(ind_valid),intensity_error_diff_avg(1,ind_valid),res) overlay(plot0,plot2) end if if (num_exp_cases.ge.3) then delete(ind_valid) ind_valid = ind(.not.ismissing(intensity_error_diff_avg(2,:))) res@xyLineColor = colors_plot(3) res@xyLineThicknessF = 8.0 plot3 = gsn_csm_xy(wks,x_axis(ind_valid),intensity_error_diff_avg(2,ind_valid),res) overlay(plot0,plot3) end if if (num_exp_cases.ge.4) then delete(ind_valid) ind_valid = ind(.not.ismissing(intensity_error_diff_avg(3,:))) res@xyLineColor = colors_plot(4) res@xyLineThicknessF = 8.0 plot4 = gsn_csm_xy(wks,x_axis(ind_valid),intensity_error_diff_avg(3,ind_valid),res) overlay(plot0,plot4) end if ; Attach a legend marker_value = (/3,4,6,9/) mres = True txres = True txres@txFontHeightF = 0.016 do k=1,dimsizes(labels2)-1 ; start at one to skip CTL label txres@txFontColor = colors_plot(k) gsn_text_ndc(wks,labels2(k),xtxt3(k),ytxt3(k),txres) if (k.ge.1) then mres@gsMarkerIndex = marker_value(k-1) mres@gsMarkerSizeF = 16.0 mres@gsMarkerThicknessF = 4.0 mres@gsMarkerColor = colors_plot(k) gsn_polymarker_ndc(wks,xtxt3(k)-0.09,ytxt3(k),mres) end if end do ; overlay 95% confidence interval polyres = True polyres@gsMarkerIndex = 1 polyres@gsMarkerSizeF = 0.04 bar_length = (/0.4,0.3,0.2,0.1/) do k=0,dimsizes(EXP_name)-1 ; loop over experiments polyres@gsLineColor = colors_plot(k+1) polyres@gsLineThicknessF = 4.0 do w=0,dimsizes(x_axis)-1 if (.not.ismissing(intensity_error_diff_avg(k,w))) then gsn_polyline(wks,plot0,(/x_axis(w),x_axis(w)/),(/1.0*interval_intensity_diff(k,w),-1.0*interval_intensity_diff(k,w)/),polyres) gsn_polyline(wks,plot0,(/x_axis(w)-bar_length(k),x_axis(w)+bar_length(k)/),(/1.0*interval_intensity_diff(k,w),1.0*interval_intensity_diff(k,w)/),polyres) ; error bar along x gsn_polyline(wks,plot0,(/x_axis(w)-bar_length(k),x_axis(w)+bar_length(k)/),(/-1.0*interval_intensity_diff(k,w),-1.0*interval_intensity_diff(k,w)/),polyres) ; error bar along x end if end do end do draw(plot0) frame(wks) delete(res) delete(polyres) ; track error max_diff = max(track_error_diff_avg) ; set bounds min_diff = min(track_error_diff_avg) range_diff = max_diff - min_diff wks = gsn_open_wks (output_type,"track_confidence_interval_fcst_hour-all-storms_" + region_file) res = True ; plot mods desired res@gsnXYBarChart = False res@tiMainString = "Track Error Difference (" + count_storms + " TCs)" + " " + region_name res@trYMaxF = max_diff + 1.2*range_diff res@trYMinF = min_diff - 1.2*range_diff res@gsnFrame = False ; don't advance frame yet res@gsnDraw = False res@gsnFrame = False res@vpXF = 0.15 res@vpYF = 0.90 res@vpWidthF = 0.7 res@vpHeightF = 0.6 res@gsnScale = True res@trXMinF = 0 res@trXMaxF = 28 x_axis = ispan(0,28,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour res@tmXBValues = ispan(0,28,2) res@tmXBLabels = ispan(0,168,12) res@tiYAxisString = "Track Error Difference (km)" res@tiXAxisString = "Forecast lead time (hr)" res@tiXAxisFontHeightF = 0.015 res@tmYROn = False res@tmYROn = False res@tmXMajorGrid = True ; implement x grid res@tmXMajorGridThicknessF = 1.0 ; 2.0 is default res@tmXMajorGridLineDashPattern = 2 res@tmYRMinorOn = False res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 8.0 res@tiMainFontHeightF = 0.021 res@xyLineColor = colors_plot(1) res@gsnYRefLine = (/0/) ; first plot delete(ind_valid) ind_valid = ind(.not.ismissing(track_error_diff_avg(0,:))) plot0 = gsn_csm_xy(wks,x_axis(ind_valid),track_error_diff_avg(0,ind_valid),res) if (num_exp_cases.ge.2) then delete(ind_valid) ind_valid = ind(.not.ismissing(track_error_diff_avg(1,:))) res@xyLineColor = colors_plot(2) res@xyLineThicknessF = 8.0 plot2 = gsn_csm_xy(wks,x_axis(ind_valid),track_error_diff_avg(1,ind_valid),res) overlay(plot0,plot2) end if if (num_exp_cases.ge.3) then delete(ind_valid) ind_valid = ind(.not.ismissing(track_error_diff_avg(2,:))) res@xyLineColor = colors_plot(3) res@xyLineThicknessF = 8.0 plot3 = gsn_csm_xy(wks,x_axis(ind_valid),track_error_diff_avg(2,ind_valid),res) overlay(plot0,plot3) end if if (num_exp_cases.ge.4) then delete(ind_valid) ind_valid = ind(.not.ismissing(track_error_diff_avg(3,:))) res@xyLineColor = colors_plot(4) res@xyLineThicknessF = 8.0 plot4 = gsn_csm_xy(wks,x_axis(ind_valid),track_error_diff_avg(3,ind_valid),res) overlay(plot0,plot4) end if ; Attach a legend marker_value = (/3,4,6,9/) mres = True txres = True txres@txFontHeightF = 0.016 do k=1,dimsizes(labels2)-1 ; start at one to skip CTL label txres@txFontColor = colors_plot(k) gsn_text_ndc(wks,labels2(k),xtxt3(k),ytxt3(k),txres) if (k.ge.1) then mres@gsMarkerIndex = marker_value(k-1) mres@gsMarkerSizeF = 16.0 mres@gsMarkerThicknessF = 4.0 mres@gsMarkerColor = colors_plot(k) gsn_polymarker_ndc(wks,xtxt3(k)-0.09,ytxt3(k),mres) end if end do ; overlay 95% confidence interval polyres = True polyres@gsMarkerIndex = 1 polyres@gsMarkerSizeF = 0.04 bar_length = (/0.4,0.3,0.2,0.1/) do k=0,dimsizes(EXP_name)-1 ; loop over experiments polyres@gsLineColor = colors_plot(k+1) polyres@gsLineThicknessF = 4.0 do w=0,dimsizes(x_axis)-1 if (.not.ismissing(track_error_diff_avg(k,w))) then gsn_polyline(wks,plot0,(/x_axis(w),x_axis(w)/),(/1.0*interval_track_diff(k,w),-1.0*interval_track_diff(k,w)/),polyres) gsn_polyline(wks,plot0,(/x_axis(w)-bar_length(k),x_axis(w)+bar_length(k)/),(/1.0*interval_track_diff(k,w),1.0*interval_track_diff(k,w)/),polyres) ; error bar along x gsn_polyline(wks,plot0,(/x_axis(w)-bar_length(k),x_axis(w)+bar_length(k)/),(/-1.0*interval_track_diff(k,w),-1.0*interval_track_diff(k,w)/),polyres) ; error bar along x end if end do end do draw(plot0) frame(wks) delete(res) delete(polyres) ; wind error max_diff = max(wind_error_diff_avg) ; set bounds min_diff = min(wind_error_diff_avg) range_diff = max_diff - min_diff wks = gsn_open_wks (output_type,"wind_confidence_interval_fcst_hour-all-storms_" + region_file) res = True ; plot mods desired res@gsnXYBarChart = False res@tiMainString = "Wind Error Difference (" + count_storms + " TCs)" + " " + region_name res@trYMaxF = max_diff + 1.2*range_diff res@trYMinF = min_diff - 1.2*range_diff res@gsnFrame = False ; don't advance frame yet res@gsnDraw = False res@gsnFrame = False res@vpXF = 0.15 res@vpYF = 0.90 res@vpWidthF = 0.7 res@vpHeightF = 0.6 res@gsnScale = True res@trXMinF = 0 res@trXMaxF = 28 x_axis = ispan(0,28,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour res@tmXBValues = ispan(0,28,2) res@tmXBLabels = ispan(0,168,12) res@tiYAxisString = "Wind Error Difference (kts)" res@tiXAxisString = "Forecast lead time (hr)" res@tiXAxisFontHeightF = 0.015 res@tmYROn = False res@tmYROn = False res@tmXMajorGrid = True ; implement x grid res@tmXMajorGridThicknessF = 1.0 ; 2.0 is default res@tmXMajorGridLineDashPattern = 2 res@tmYRMinorOn = False res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 8.0 res@tiMainFontHeightF = 0.021 res@xyLineColor = colors_plot(1) res@gsnYRefLine = (/0/) ; first plot delete(ind_valid) ind_valid = ind(.not.ismissing(wind_error_diff_avg(0,:))) plot0 = gsn_csm_xy(wks,x_axis(ind_valid),wind_error_diff_avg(0,ind_valid),res) if (num_exp_cases.ge.2) then delete(ind_valid) ind_valid = ind(.not.ismissing(wind_error_diff_avg(1,:))) res@xyLineColor = colors_plot(2) res@xyLineThicknessF = 8.0 plot2 = gsn_csm_xy(wks,x_axis(ind_valid),wind_error_diff_avg(1,ind_valid),res) overlay(plot0,plot2) end if if (num_exp_cases.ge.3) then delete(ind_valid) ind_valid = ind(.not.ismissing(wind_error_diff_avg(2,:))) res@xyLineColor = colors_plot(3) res@xyLineThicknessF = 8.0 plot3 = gsn_csm_xy(wks,x_axis(ind_valid),wind_error_diff_avg(2,ind_valid),res) overlay(plot0,plot3) end if if (num_exp_cases.ge.4) then delete(ind_valid) ind_valid = ind(.not.ismissing(wind_error_diff_avg(3,:))) res@xyLineColor = colors_plot(4) res@xyLineThicknessF = 8.0 plot4 = gsn_csm_xy(wks,x_axis(ind_valid),wind_error_diff_avg(3,ind_valid),res) overlay(plot0,plot4) end if ; Attach a legend marker_value = (/3,4,6,9/) mres = True txres = True txres@txFontHeightF = 0.016 do k=1,dimsizes(labels2)-1 ; start at one to skip CTL label txres@txFontColor = colors_plot(k) gsn_text_ndc(wks,labels2(k),xtxt3(k),ytxt3(k),txres) if (k.ge.1) then mres@gsMarkerIndex = marker_value(k-1) mres@gsMarkerSizeF = 16.0 mres@gsMarkerThicknessF = 4.0 mres@gsMarkerColor = colors_plot(k) gsn_polymarker_ndc(wks,xtxt3(k)-0.09,ytxt3(k),mres) end if end do ; overlay 95% confidence interval polyres = True polyres@gsMarkerIndex = 1 polyres@gsMarkerSizeF = 0.04 bar_length = (/0.4,0.3,0.2,0.1/) do k=0,dimsizes(EXP_name)-1 ; loop over experiments polyres@gsLineColor = colors_plot(k+1) polyres@gsLineThicknessF = 4.0 do w=0,dimsizes(x_axis)-1 if (.not.ismissing(wind_error_diff_avg(k,w))) then gsn_polyline(wks,plot0,(/x_axis(w),x_axis(w)/),(/1.0*interval_wind_diff(k,w),-1.0*interval_wind_diff(k,w)/),polyres) gsn_polyline(wks,plot0,(/x_axis(w)-bar_length(k),x_axis(w)+bar_length(k)/),(/1.0*interval_wind_diff(k,w),1.0*interval_wind_diff(k,w)/),polyres) ; error bar along x gsn_polyline(wks,plot0,(/x_axis(w)-bar_length(k),x_axis(w)+bar_length(k)/),(/-1.0*interval_wind_diff(k,w),-1.0*interval_wind_diff(k,w)/),polyres) ; error bar along x end if end do end do draw(plot0) frame(wks) delete(res) delete(polyres) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; delete([/track_temp,max_track,min_track,max_int,min_int,int_temp,max_fcst,min_fcst,x_axis,check_exp,track_temp,int_temp,wind_temp/]) check_exp = new((/dimsizes(EXP_name)/),"float") count_exp1 = new((/dimsizes(EXP_name)/),"integer") ;;;;;;;;;;;;;;; ;;;; plot average track and intensity errors as function of forecast hour over all storms individually ;;;;;;;;;;;;;;; track_error_std_st_fh = new((/dimsizes(EXP_name),sizes_g5nr,29/),"float") intensity_error_std_st_fh = track_error_std_st_fh wind_error_std_st_fh = track_error_std_st_fh interval_track_st_fh = track_error_std_st_fh interval_intensity_st_fh = track_error_std_st_fh interval_wind_st_fh = track_error_std_st_fh track_error_avg_st_fh = track_error_std_st_fh intensity_error_avg_st_fh = track_error_std_st_fh wind_error_avg_st_fh = track_error_std_st_fh max_ctl_length_st = dim_num_n(.not.ismissing(track_error_ctl_st_fh),(/1/)) ; max length for each paired experiment and storm ; 95 % confidence inteval for individual storms do a=0,dimsizes(EXP_name)-1 track_error_std_st_fh(a,:,:) = dim_stddev_n(track_error_exp_conform(a,:,:,:)-track_error_ctl_conform(a,:,:,:),1) ; forecast hour only intensity_error_std_st_fh(a,:,:) = dim_stddev_n(intensity_error_exp_conform(a,:,:,:)-intensity_error_ctl_conform(a,:,:,:),1) ; forecast hour only, wind_error_std_st_fh(a,:,:) = dim_stddev_n(wind_error_exp_conform(a,:,:,:)-wind_error_ctl_conform(a,:,:,:),1) ; forecast hour only track_error_avg_st_fh(a,:,:) = dim_avg_n(track_error_exp_conform(a,:,:,:)-track_error_ctl_conform(a,:,:,:),1) ; forecast hour only intensity_error_avg_st_fh(a,:,:) = dim_avg_n(intensity_error_exp_conform(a,:,:,:)-intensity_error_ctl_conform(a,:,:,:),1) ; forecast hour only wind_error_avg_st_fh(a,:,:) = dim_avg_n(wind_error_exp_conform(a,:,:,:)-wind_error_ctl_conform(a,:,:,:),1) ; forecast hour only, end do ; end do here over experiment names print("Making average results as function of forecast hour for each storm...") do i=0,dimsizes(g5nr_name)-1 ; loop over all storms ; check various experiments do a=0,dimsizes(EXP_name)-1 if (all(ismissing(track_error_exp_st_fh(i,:,a)))) then check_exp(a) = 1.0 else check_exp(a) = 0.0 end if end do location_max_st_temp = ind(max_ctl_length_st(i,:).eq.max(max_ctl_length_st(i,:))) ; max ctl length for this storm location_max_st = location_max_st_temp(0) ; delete this at end ; check to see if there is missing data, meaning no data for this storm if (.not.all(ismissing(track_error_ctl_st_fh(i,:,location_max_st))) .and. all(check_exp.eq.0.0) ) then ; then make track and intensity error plots if (num(.not.ismissing(track_error_ctl_st_fh(i,:,location_max_st))).gt.1) then ; need 2 elements to plot a line ; then make plots print("Making plot for: " + g5nr_name(i)) track_temp = array_append_record(track_error_ctl_st_fh(i,:,location_max_st),ndtooned(track_error_exp_st_fh(i,:,:)),0) max_track = max(track_temp) min_track = min(track_temp) int_temp = array_append_record(intensity_error_ctl_st_fh(i,:,location_max_st),ndtooned(intensity_error_exp_st_fh(i,:,:)),0) max_int = max(int_temp) min_int = min(int_temp) wind_temp = array_append_record(wind_error_ctl_st_fh(i,:,location_max_st),ndtooned(wind_error_exp_st_fh(i,:,:)),0) max_wind = max(wind_temp) min_wind = min(wind_temp) ; determine number of samples do a=0,dimsizes(EXP_name)-1 do k=0,28 ; loop over forecast hours (output every 6 hours) fcst_hr_array_exp2(a,k) = num(.not.ismissing(track_error_exp_conform(a,i,:,k))) fcst_hr_array_ctl2(a,k) = num(.not.ismissing(track_error_ctl_conform(a,i,:,k))) end do end do max_fcst = max(fcst_hr_array_ctl2) min_fcst = min(fcst_hr_array_ctl2) ; count number of non-missing values, conform each to the smaller array do a=0,dimsizes(EXP_name)-1 count_exp1(a) = num(.not.ismissing(track_error_exp_st_fh(i,:,a))) end do min_countexp1 = min(count_exp1) count_ctl = num(.not.ismissing(track_error_ctl_st_fh(i,:,location_max_st))) if (count_ctl.le.min_countexp1) then x_axis_max = tointeger(count_ctl) - 1 ; subtract one based on array indexing end if if (min_countexp1.lt.count_ctl) then x_axis_max = tointeger(min_countexp1) - 1 ; subtract one based on array indexing end if do a=0,dimsizes(EXP_name)-1 do y=0,28 ; loop over samples to divide by smallest of the two samples if (fcst_hr_array_ctl2(a,y).le.fcst_hr_array_exp2(a,y)) then number_sample(a,y) = fcst_hr_array_ctl2(a,y) if (number_sample(a,y).lt.20) then const_z(y) = 2.228 end if if (number_sample(a,y).ge.20 .and. number_sample(a,y).lt.40) then const_z(y) = 2.042 end if if (number_sample(a,y).ge.40 .and. number_sample(a,y).lt.80) then const_z(y) = 2.000 end if if (number_sample(a,y).ge.80) then const_z(y) = 1.96 end if else number_sample(a,y) = fcst_hr_array_exp2(a,y) if (number_sample(a,y).lt.20) then const_z(y) = 2.228 end if if (number_sample(a,y).ge.20 .and. number_sample(a,y).lt.40) then const_z(y) = 2.042 end if if (number_sample(a,y).ge.40 .and. number_sample(a,y).lt.80) then const_z(y) = 2.000 end if if (number_sample(a,y).ge.80) then const_z(y) = 1.96 end if end if end do end do do a=0,dimsizes(EXP_name)-1 number_valid = ind(number_sample(a,:).ne.0) if (dimsizes(number_valid).ge.1) then interval_track_st_fh(a,i,number_valid) = const_z(number_valid)*track_error_std_st_fh(a,i,number_valid) / sqrt(number_sample(a,number_valid)) ; % 95 interval interval_intensity_st_fh(a,i,number_valid) = const_z(number_valid)*intensity_error_std_st_fh(a,i,number_valid) / sqrt(number_sample(a,number_valid)) ; 95% interval interval_wind_st_fh(a,i,number_valid) = const_z(number_valid)*wind_error_std_st_fh(a,i,number_valid) / sqrt(number_sample(a,number_valid)) ; 95% interval end if delete(number_valid) end do ;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; plot errors as a function of forecast lead time ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; wks = gsn_open_wks (output_type,"track_error_fcst_hour-avg-" + g5nr_name(i)) res = True ; plot mods desired res@tiMainString = "Track Error (" + g5nr_name(i) + ")" ; plot mods desired res@gsnXYBarChart = False res@trYMaxF = max_track + 100.0 res@trYMinF = 0 res@gsnFrame = False ; don't advance frame yet res@gsnDraw = False res@gsnFrame = False res@vpXF = 0.15 res@vpYF = 0.90 res@vpWidthF = 0.7 res@vpHeightF = 0.6 res@gsnScale = True res@trXMinF = 0 res@trXMaxF = 28 x_axis = ispan(0,28,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour res@tmXBValues = ispan(0,28,2) res@tmXBLabels = ispan(0,168,12) ; x_axis = ispan(0,x_axis_max,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour ; res@tmXBValues = ispan(0,x_axis_max,2) ; res@tmXBLabels = fcst_hr(0:x_axis_max:2) res@tiYAxisString = "Track Error (km)" res@tiXAxisString = "Forecast lead time (hr)" res@tiXAxisFontHeightF = 0.015 res@tmYROn = False res@tmYROn = False res@tmYRMinorOn = False res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 8.0 res@tmXMajorGrid = True ; implement x grid res@tmXMajorGridThicknessF = 1.0 ; 2.0 is default res@tmXMajorGridLineDashPattern = 2 res@tiMainFontHeightF = 0.021 res@xyLineColor = colors_plot(0) plot0 = gsn_csm_xy(wks,x_axis(0:count_ctl-1),track_error_ctl_st_fh(i,0:count_ctl-1,location_max_st),res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 8.0 plot1 = gsn_csm_xy(wks,x_axis(0:count_exp1(0)-1),track_error_exp_st_fh(i,0:count_exp1(0)-1,0),res) if (num_exp_cases.ge.2) then res@xyLineColor = colors_plot(2) res@xyLineThicknessF = 8.0 plot2 = gsn_csm_xy(wks,x_axis(0:count_exp1(1)-1),track_error_exp_st_fh(i,0:count_exp1(1)-1,1),res) overlay(plot0,plot2) end if if (num_exp_cases.ge.3) then res@xyLineColor = colors_plot(3) res@xyLineThicknessF = 8.0 plot3 = gsn_csm_xy(wks,x_axis(0:count_exp1(2)-1),track_error_exp_st_fh(i,0:count_exp1(2)-1,2),res) overlay(plot0,plot3) end if if (num_exp_cases.ge.4) then res@xyLineColor = colors_plot(4) res@xyLineThicknessF = 8.0 plot4 = gsn_csm_xy(wks,x_axis(0:count_exp1(3)-1),track_error_exp_st_fh(i,0:count_exp1(3)-1,3),res) overlay(plot0,plot4) end if bres = True bres@tmXBMode = "Explicit" bres@gsnFrame = False bres@gsnDraw = False bres@gsnXYBarChart = True bres@gsnXYBarChartColors = "grey78" bres@tmYROn = True bres@tmYRMinorOn = True bres@vpXF = 0.15 bres@vpYF = 0.90 bres@vpWidthF = 0.7 bres@vpHeightF = 0.6 bres@tmYLOn = False bres@trYMinF = 0 bres@tiMainString = "" bres@trYMaxF = max_fcst + 4.0 ; number of samples in each forecast hour bres@tmYRLabelsOn = True bres@tiYAxisFontHeightF = 0.015 bres@tmYRLabelFontHeightF = 0.015 bres@gsnXYBarChartOutlineThicknessF = 1.5 bres@tmYROn = True bres@gsnXYBarChartBarWidth = 0.5 ; change bar widths 0.75 originally bres@tmXBOn = False bres@trXMinF = 0 ; adds space on either end bres@trXMaxF = x_axis_max ; of the 1st and last bars bres@tiYAxisString = "# fcsts" bres@tiYAxisSide = "Right" bres@tiYAxisAngleF = 270 bres@tmXTOn = False ;plot5 = gsn_csm_xy(wks,x_axis,fcst_hr_array_ctl2(location_max_st,0:x_axis_max),bres) bres@gsnXYBarChartColors = "grey78" ;plot6 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp2(0,0:x_axis_max),bres) ;draw(plot5) ;draw(plot6) if (num_exp_cases.ge.2) then bres@gsnXYBarChartColors = "grey78" ;plot7 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp2(1,0:x_axis_max),bres) ;draw(plot7) end if if (num_exp_cases.ge.3) then bres@gsnXYBarChartColors = "grey78" ;plot8 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp2(2,0:x_axis_max),bres) ;draw(plot8) end if if (num_exp_cases.ge.4) then bres@gsnXYBarChartColors = "grey78" ;plot9 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp2(3,0:x_axis_max),bres) ;draw(plot9) end if overlay(plot0,plot1) ; Attach a legend marker_value = (/3,4,6,9/) mres = True txres = True txres@txFontHeightF = 0.016 do k=0,dimsizes(labels2)-1 txres@txFontColor = colors_plot(k) gsn_text_ndc(wks,labels2(k),xtxt3(k),ytxt3(k),txres) if (k.ge.1) then mres@gsMarkerIndex = marker_value(k-1) mres@gsMarkerSizeF = 16.0 mres@gsMarkerThicknessF = 4.0 mres@gsMarkerColor = colors_plot(k) gsn_polymarker_ndc(wks,xtxt3(k)-0.07,ytxt3(k),mres) end if end do ; attach significance do k=0,dimsizes(EXP_name)-1 do w=0,x_axis_max ; loop over forecast hours mres@gsMarkerIndex = marker_value(k) mres@gsMarkerSizeF = 16.0 mres@gsMarkerThicknessF = 4.0 ;mres@gsMarkerColor = colors_plot(k+1) if (.not.ismissing(track_error_avg_st_fh(k,i,w)) .and. interval_track_st_fh(k,i,w).gt.0) then if (abs(track_error_avg_st_fh(k,i,w)).gt.interval_track_st_fh(k,i,w)) then ; 95% confidence interval if (track_error_avg_st_fh(k,i,w).lt.0) then ; then exp is better mres@gsMarkerColor = colors_plot(k+1) else ; ctl is better mres@gsMarkerColor = colors_plot(0) end if gsn_polymarker(wks,plot0,x_axis(w),10.0,mres) end if end if end do end do ;;;;; resbot = True resbot@vpXF = 0.15 ; The left side of the box location resbot@vpYF = 0.20 ; The top side of the plot box location resbot@vpWidthF = 0.70 ; The Width of the plot box resbot@vpHeightF = 0.05 ; The height of the plot box resbot@trYMaxF = 3 ; random value does not need to be changed resbot@trXMinF = 0 resbot@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour resbot@tiYAxisString = "# fcsts" resbot@trXMaxF = 28 resbot@tmXBValues = ispan(0,28,2) resbot@tmXBLabels = ispan(0,168,12) resbot@tiXAxisString = "Forecast lead time (hr)" resbot@tiXAxisFontHeightF = 0.015 resbot@tmYROn = False resbot@tmYROn = False resbot@tmYRMinorOn = False resbot@tmYLBorderOn = False resbot@tmYLMinorPerMajor = 1 resbot@tmYLOn = False resbot@tmYRBorderOn = False resbot@tmYROn = False resbot@tiYAxisFontHeightF = 0.015 resbot@tmXBLabelFontHeightF = 0.015 resbot@tmYLLabelFontHeightF = 0.015 resbot@xyLineThicknessF = 0.05 resbot@tiMainFontHeightF = 0.021 resbot@tmXTOn = False ; turn off the top tickmarks resbot@tmYLOn = False resbot@tmXBMinorOn = False ; No minor tick marks. resbot@gsnDraw = False ; Don't draw individual plot. resbot@gsnFrame = False ; Don't advance frame. ;resbot@gsnYRefLine = 0.0 ; create a reference line ;resbot@gsnAboveYRefLineColor = "green" ; above ref line fill green resbot@gsnXYBarChart = False ; turn on bar chart plot_blank = gsn_csm_xy(wks,(/1,2,3/),(/0.0,0.0,0.0/),resbot) ; plot false values to create plot, then overlay text values on plot ; attache the number of paired cases over the plot y_value_text = (/1.5,1.5,1.5/) do k=0,dimsizes(EXP_name)-1 do w=0,28,2 ; loop over forecast hours, skip first hour for simplicity of viewing image qres = True qres@txFontHeightF = 0.017 qres@txJust = "CenterCenter" qres@txFontColor = colors_plot(k+1) ; color coded by experiment text = gsn_add_text(wks,plot_blank,str_squeeze(tostring(tointeger(fcst_hr_array_exp2(k,w)))),x_axis(w),y_value_text(k),qres) end do end do ;;;;;;;;;;;;; draw(plot0) draw(plot_blank) frame(wks) delete(res) delete(bres) delete(resbot) ;;;;;;;;;;;;;;;;;;; wks = gsn_open_wks (output_type,"intensity_error_fcst_hour-avg-" + g5nr_name(i)) res = True ; plot mods desired res@tiMainString = "Intensity Error (" + g5nr_name(i) + ")" ; plot mods desired res@gsnXYBarChart = False res@trYMaxF = max_int + 20.0 res@trYMinF = 0 res@gsnFrame = False ; don't advance frame yet res@gsnDraw = False res@gsnFrame = False res@vpXF = 0.15 res@vpYF = 0.90 res@vpWidthF = 0.7 res@vpHeightF = 0.6 res@gsnScale = True res@trXMinF = 0 res@trXMaxF = 28 x_axis = ispan(0,28,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour res@tmXBValues = ispan(0,28,2) res@tmXBLabels = ispan(0,168,12) ; x_axis = ispan(0,x_axis_max,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour ; res@tmXBValues = ispan(0,x_axis_max,2) ; res@tmXBLabels = fcst_hr(0:x_axis_max:2) res@tiYAxisString = "Intensity Error (hPa)" res@tiXAxisString = "Forecast lead time (hr)" res@tiXAxisFontHeightF = 0.015 res@tmYROn = False res@tmYROn = False res@tmYRMinorOn = False res@tmXMajorGrid = True ; implement x grid res@tmXMajorGridThicknessF = 1.0 ; 2.0 is default res@tmXMajorGridLineDashPattern = 2 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 8.0 res@tiMainFontHeightF = 0.021 res@xyLineColor = colors_plot(0) plot0 = gsn_csm_xy(wks,x_axis(0:count_ctl-1),intensity_error_ctl_st_fh(i,0:count_ctl-1,location_max_st),res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 8.0 plot1 = gsn_csm_xy(wks,x_axis(0:count_exp1(0)-1),intensity_error_exp_st_fh(i,0:count_exp1(0)-1,0),res) if (num_exp_cases.ge.2) then res@xyLineColor = colors_plot(2) res@xyLineThicknessF = 8.0 plot2 = gsn_csm_xy(wks,x_axis(0:count_exp1(1)-1),intensity_error_exp_st_fh(i,0:count_exp1(1)-1,1),res) overlay(plot0,plot2) end if if (num_exp_cases.ge.3) then res@xyLineColor = colors_plot(3) res@xyLineThicknessF = 8.0 plot3 = gsn_csm_xy(wks,x_axis(0:count_exp1(2)-1),intensity_error_exp_st_fh(i,0:count_exp1(2)-1,2),res) overlay(plot0,plot3) end if if (num_exp_cases.ge.4) then res@xyLineColor = colors_plot(4) res@xyLineThicknessF = 8.0 plot4 = gsn_csm_xy(wks,x_axis(0:count_exp1(3)-1),intensity_error_exp_st_fh(i,0:count_exp1(3)-1,3),res) overlay(plot0,plot4) end if bres = True bres@tmXBMode = "Explicit" bres@gsnFrame = False bres@gsnDraw = False bres@gsnXYBarChart = True bres@gsnXYBarChartColors = "grey78" bres@tmYROn = True bres@tmYRMinorOn = True bres@tmYLOn = False bres@vpXF = 0.15 bres@vpYF = 0.90 bres@vpWidthF = 0.7 bres@vpHeightF = 0.6 bres@trYMinF = 0 bres@tiMainString = "" bres@trYMaxF = max_fcst + 4.0 ; number of samples in each forecast hour bres@tmYRLabelsOn = True bres@tiYAxisFontHeightF = 0.015 bres@tmYRLabelFontHeightF = 0.015 bres@gsnXYBarChartOutlineThicknessF = 1.5 bres@tmYROn = True bres@gsnXYBarChartBarWidth = 0.5 ; change bar widths 0.75 originally bres@tmXBOn = False bres@trXMinF = 0 ; adds space on either end bres@trXMaxF = x_axis_max ; of the 1st and last bars bres@tiYAxisString = "# fcsts" bres@tiYAxisSide = "Right" bres@tiYAxisAngleF = 270 bres@tmXTOn = False ;plot5 = gsn_csm_xy(wks,x_axis,fcst_hr_array_ctl2(location_max_st,0:x_axis_max),bres) bres@gsnXYBarChartColors = "grey78" ;plot6 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp2(0,0:x_axis_max),bres) ;draw(plot5) ;draw(plot6) if (num_exp_cases.ge.2) then bres@gsnXYBarChartColors = "grey78" ;plot7 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp2(1,0:x_axis_max),bres) ;draw(plot7) end if if (num_exp_cases.ge.3) then bres@gsnXYBarChartColors = "grey78" ;plot8 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp2(2,0:x_axis_max),bres) ;draw(plot8) end if if (num_exp_cases.ge.4) then bres@gsnXYBarChartColors = "grey78" ;plot9 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp2(3,0:x_axis_max),bres) ;draw(plot9) end if overlay(plot0,plot1) ; Attach a legend marker_value = (/3,4,6,9/) mres = True txres = True txres@txFontHeightF = 0.016 do k=0,dimsizes(labels2)-1 txres@txFontColor = colors_plot(k) gsn_text_ndc(wks,labels2(k),xtxt3(k),ytxt3(k),txres) if (k.ge.1) then mres@gsMarkerIndex = marker_value(k-1) mres@gsMarkerSizeF = 16.0 mres@gsMarkerThicknessF = 4.0 mres@gsMarkerColor = colors_plot(k) gsn_polymarker_ndc(wks,xtxt3(k)-0.07,ytxt3(k),mres) end if end do ; attach significance do k=0,dimsizes(EXP_name)-1 do w=0,x_axis_max ; loop over forecast hours mres@gsMarkerIndex = marker_value(k) mres@gsMarkerSizeF = 16.0 mres@gsMarkerThicknessF = 4.0 ;mres@gsMarkerColor = colors_plot(k+1) if (.not.ismissing(intensity_error_avg_st_fh(k,i,w)) .and. interval_intensity_st_fh(k,i,w).gt.0) then if (abs(intensity_error_avg_st_fh(k,i,w)).gt.interval_intensity_st_fh(k,i,w)) then ; 95% confidence interval if (intensity_error_avg_st_fh(k,i,w).lt.0) then ; then exp is better mres@gsMarkerColor = colors_plot(k+1) else ; ctl is better mres@gsMarkerColor = colors_plot(0) end if gsn_polymarker(wks,plot0,x_axis(w),2.0,mres) end if end if end do end do ;;;;; resbot = True resbot@vpXF = 0.15 ; The left side of the box location resbot@vpYF = 0.20 ; The top side of the plot box location resbot@vpWidthF = 0.70 ; The Width of the plot box resbot@vpHeightF = 0.05 ; The height of the plot box resbot@trYMaxF = 3 ; random value does not need to be changed resbot@trXMinF = 0 resbot@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour resbot@tiYAxisString = "# fcsts" resbot@tiXAxisString = "Forecast lead time (hr)" resbot@tiXAxisFontHeightF = 0.015 resbot@trXMaxF = 28 resbot@tmXBValues = ispan(0,28,2) resbot@tmXBLabels = ispan(0,168,12) resbot@tmYROn = False resbot@tmYROn = False resbot@tmYLBorderOn = False resbot@tmYLMinorPerMajor = 1 resbot@tmYLOn = False resbot@tmYRBorderOn = False resbot@tmYROn = False resbot@tmYRMinorOn = False resbot@tiYAxisFontHeightF = 0.015 resbot@tmXBLabelFontHeightF = 0.015 resbot@tmYLLabelFontHeightF = 0.015 resbot@xyLineThicknessF = 0.05 resbot@tiMainFontHeightF = 0.021 resbot@tmXTOn = False ; turn off the top tickmarks resbot@tmYLOn = False resbot@tmXBMinorOn = False ; No minor tick marks. resbot@gsnDraw = False ; Don't draw individual plot. resbot@gsnFrame = False ; Don't advance frame. ;resbot@gsnYRefLine = 0.0 ; create a reference line ;resbot@gsnAboveYRefLineColor = "green" ; above ref line fill green resbot@gsnXYBarChart = False ; turn on bar chart plot_blank = gsn_csm_xy(wks,(/1,2,3/),(/0.0,0.0,0.0/),resbot) ; plot false values to create plot, then overlay text values on plot ; attache the number of paired cases over the plot y_value_text = (/1.5,1.5,1.5/) do k=0,dimsizes(EXP_name)-1 do w=0,28,2 ; loop over forecast hours, skip first hour for simplicity of viewing image qres = True qres@txFontHeightF = 0.017 qres@txJust = "CenterCenter" qres@txFontColor = colors_plot(k+1) ; color coded by experiment text = gsn_add_text(wks,plot_blank,str_squeeze(tostring(tointeger(fcst_hr_array_exp2(k,w)))),x_axis(w),y_value_text(k),qres) end do end do ;;;;;;;;;;;;; draw(plot0) draw(plot_blank) frame(wks) delete(res) delete(bres) delete(resbot) ; wind wks = gsn_open_wks (output_type,"wind_error_fcst_hour-avg-" + g5nr_name(i)) res = True ; plot mods desired res@tiMainString = "Wind Error (" + g5nr_name(i) + ")" ; plot mods desired res@gsnXYBarChart = False res@trYMaxF = max_wind + 20.0 res@trYMinF = 0 res@gsnFrame = False ; don't advance frame yet res@gsnDraw = False res@gsnFrame = False res@vpXF = 0.15 res@vpYF = 0.90 res@vpWidthF = 0.7 res@vpHeightF = 0.6 res@gsnScale = True res@trXMinF = 0 res@trXMaxF = 28 x_axis = ispan(0,28,1) ;x_axis = ispan(0,x_axis_max,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour res@tmXBValues = ispan(0,28,2) res@tmXBLabels = ispan(0,168,12) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour ;res@tmXBValues = ispan(0,x_axis_max,2) ;res@tmXBLabels = fcst_hr(0:x_axis_max:2) res@tiYAxisString = "Wind Error (kts)" res@tiXAxisString = "Forecast lead time (hr)" res@tiXAxisFontHeightF = 0.015 res@tmYROn = False res@tmYROn = False res@tmYRMinorOn = False res@tmXMajorGrid = True ; implement x grid res@tmXMajorGridThicknessF = 1.0 ; 2.0 is default res@tmXMajorGridLineDashPattern = 2 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 8.0 res@tiMainFontHeightF = 0.021 res@xyLineColor = colors_plot(0) plot0 = gsn_csm_xy(wks,x_axis(0:count_ctl-1),wind_error_ctl_st_fh(i,0:count_ctl-1,location_max_st),res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 8.0 plot1 = gsn_csm_xy(wks,x_axis(0:count_exp1(0)-1),wind_error_exp_st_fh(i,0:count_exp1(0)-1,0),res) if (num_exp_cases.ge.2) then res@xyLineColor = colors_plot(2) res@xyLineThicknessF = 8.0 plot2 = gsn_csm_xy(wks,x_axis(0:count_exp1(1)-1),wind_error_exp_st_fh(i,0:count_exp1(1)-1,1),res) overlay(plot0,plot2) end if if (num_exp_cases.ge.3) then res@xyLineColor = colors_plot(3) res@xyLineThicknessF = 8.0 plot3 = gsn_csm_xy(wks,x_axis(0:count_exp1(2)-1),wind_error_exp_st_fh(i,0:count_exp1(2)-1,2),res) overlay(plot0,plot3) end if if (num_exp_cases.ge.4) then res@xyLineColor = colors_plot(4) res@xyLineThicknessF = 8.0 plot4 = gsn_csm_xy(wks,x_axis(0:count_exp1(3)-1),wind_error_exp_st_fh(i,0:count_exp1(3)-1,3),res) overlay(plot0,plot4) end if bres = True bres@tmXBMode = "Explicit" bres@gsnFrame = False bres@gsnDraw = False bres@gsnXYBarChart = True bres@gsnXYBarChartColors = "grey78" bres@tmYROn = True bres@tmYRMinorOn = True bres@vpXF = 0.15 bres@vpYF = 0.90 bres@vpWidthF = 0.7 bres@vpHeightF = 0.6 bres@tmYLOn = False bres@trYMinF = 0 bres@tiMainString = "" bres@trYMaxF = max_fcst + 4.0 ; number of samples in each forecast hour bres@tmYRLabelsOn = True bres@tiYAxisFontHeightF = 0.015 bres@tmYRLabelFontHeightF = 0.015 bres@gsnXYBarChartOutlineThicknessF = 1.5 bres@tmYROn = True bres@gsnXYBarChartBarWidth = 0.5 ; change bar widths 0.75 originally bres@tmXBOn = False bres@trXMinF = 0 ; adds space on either end bres@trXMaxF = x_axis_max ; of the 1st and last bars bres@tiYAxisString = "# fcsts" bres@tiYAxisSide = "Right" bres@tiYAxisAngleF = 270 bres@tmXTOn = False ;plot5 = gsn_csm_xy(wks,x_axis,fcst_hr_array_ctl2(location_max_st,0:x_axis_max),bres) bres@gsnXYBarChartColors = "grey78" ;plot6 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp2(0,0:x_axis_max),bres) ;draw(plot5) ;draw(plot6) if (num_exp_cases.ge.2) then bres@gsnXYBarChartColors = "grey78" ;plot7 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp2(1,0:x_axis_max),bres) ;draw(plot7) end if if (num_exp_cases.ge.3) then bres@gsnXYBarChartColors = "grey78" ;plot8 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp2(2,0:x_axis_max),bres) ;draw(plot8) end if if (num_exp_cases.ge.4) then bres@gsnXYBarChartColors = "grey78" ;plot9 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp2(3,0:x_axis_max),bres) ;draw(plot9) end if overlay(plot0,plot1) ; Attach a legend marker_value = (/3,4,6,9/) mres = True txres = True txres@txFontHeightF = 0.016 do k=0,dimsizes(labels2)-1 txres@txFontColor = colors_plot(k) gsn_text_ndc(wks,labels2(k),xtxt3(k),ytxt3(k),txres) if (k.ge.1) then mres@gsMarkerIndex = marker_value(k-1) mres@gsMarkerSizeF = 16.0 mres@gsMarkerThicknessF = 4.0 mres@gsMarkerColor = colors_plot(k) gsn_polymarker_ndc(wks,xtxt3(k)-0.07,ytxt3(k),mres) end if end do ; attach significance do k=0,dimsizes(EXP_name)-1 do w=0,x_axis_max ; loop over forecast hours mres@gsMarkerIndex = marker_value(k) mres@gsMarkerSizeF = 16.0 mres@gsMarkerThicknessF = 4.0 ;mres@gsMarkerColor = colors_plot(k+1) if (.not.ismissing(wind_error_avg_st_fh(k,i,w)) .and. interval_wind_st_fh(k,i,w).gt.0) then if (abs(wind_error_avg_st_fh(k,i,w)).gt.interval_wind_st_fh(k,i,w)) then ; 95% confidence interval if (wind_error_avg_st_fh(k,i,w).lt.0) then ; then exp is better mres@gsMarkerColor = colors_plot(k+1) else ; ctl is better mres@gsMarkerColor = colors_plot(0) end if gsn_polymarker(wks,plot0,x_axis(w),2.0,mres) end if end if end do end do ;;;;; resbot = True resbot@vpXF = 0.15 ; The left side of the box location resbot@vpYF = 0.20 ; The top side of the plot box location resbot@vpWidthF = 0.70 ; The Width of the plot box resbot@vpHeightF = 0.05 ; The height of the plot box resbot@trYMaxF = 3 ; random value does not need to be changed resbot@trXMinF = 0 resbot@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour resbot@tiYAxisString = "# fcsts" resbot@tiXAxisString = "Forecast lead time (hr)" resbot@tiXAxisFontHeightF = 0.015 resbot@trXMaxF = 28 resbot@tmXBValues = ispan(0,28,2) resbot@tmXBLabels = ispan(0,168,12) resbot@tmYROn = False resbot@tmYROn = False resbot@tmYLBorderOn = False resbot@tmYLMinorPerMajor = 1 resbot@tmYLOn = False resbot@tmYRBorderOn = False resbot@tmYROn = False resbot@tmYRMinorOn = False resbot@tiYAxisFontHeightF = 0.015 resbot@tmXBLabelFontHeightF = 0.015 resbot@tmYLLabelFontHeightF = 0.015 resbot@xyLineThicknessF = 0.05 resbot@tiMainFontHeightF = 0.021 resbot@tmXTOn = False ; turn off the top tickmarks resbot@tmYLOn = False resbot@tmXBMinorOn = False ; No minor tick marks. resbot@gsnDraw = False ; Don't draw individual plot. resbot@gsnFrame = False ; Don't advance frame. ;resbot@gsnYRefLine = 0.0 ; create a reference line ;resbot@gsnAboveYRefLineColor = "green" ; above ref line fill green resbot@gsnXYBarChart = False ; turn on bar chart plot_blank = gsn_csm_xy(wks,(/1,2,3/),(/0.0,0.0,0.0/),resbot) ; plot false values to create plot, then overlay text values on plot ; attache the number of paired cases over the plot y_value_text = (/1.5,1.5,1.5/) do k=0,dimsizes(EXP_name)-1 do w=0,28,2 ; loop over forecast hours, skip first hour for simplicity of viewing image qres = True qres@txFontHeightF = 0.017 qres@txJust = "CenterCenter" qres@txFontColor = colors_plot(k+1) ; color coded by experiment text = gsn_add_text(wks,plot_blank,str_squeeze(tostring(tointeger(fcst_hr_array_exp2(k,w)))),x_axis(w),y_value_text(k),qres) end do end do ;;;;;;;;;;;;; draw(plot0) draw(plot_blank) frame(wks) delete(res) delete(bres) delete(resbot) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; 95% CONFIDENCE INTERVAL DIFFERNCE PLOTS FOR TRACK, WIND, AND INTENSITY ; track error max_diff = max(track_error_avg_st_fh(:,i,:)) ; set bounds min_diff = min(track_error_avg_st_fh(:,i,:)) range_diff = max_diff - min_diff wks = gsn_open_wks (output_type,"track_confidence_interval_fcst_hour-avg-" + g5nr_name(i)) res = True ; plot mods desired res@gsnXYBarChart = False res@tiMainString = "Track Error Difference (" + g5nr_name(i) + ")" res@trYMaxF = max_diff + 1.2*range_diff res@trYMinF = min_diff - 1.2*range_diff res@gsnFrame = False ; don't advance frame yet res@gsnDraw = False res@gsnFrame = False res@vpXF = 0.15 res@vpYF = 0.90 res@vpWidthF = 0.7 res@vpHeightF = 0.6 res@gsnScale = True res@trXMinF = 0 res@trXMaxF = 28 x_axis = ispan(0,28,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour res@tmXBValues = ispan(0,28,2) res@tmXBLabels = ispan(0,168,12) res@tiYAxisString = "Track Error Difference (km)" res@tiXAxisString = "Forecast lead time (hr)" res@tiXAxisFontHeightF = 0.015 res@tmYROn = False res@tmYROn = False res@tmXMajorGrid = True ; implement x grid res@tmXMajorGridThicknessF = 1.0 ; 2.0 is default res@tmXMajorGridLineDashPattern = 2 res@tmYRMinorOn = False res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 8.0 res@tiMainFontHeightF = 0.021 res@xyLineColor = colors_plot(1) res@gsnYRefLine = (/0/) ; first plot plot0 = gsn_csm_xy(wks,x_axis(0:count_exp1(0)-1),track_error_avg_st_fh(0,i,0:count_exp1(0)-1),res) if (num_exp_cases.ge.2) then res@xyLineColor = colors_plot(2) res@xyLineThicknessF = 8.0 plot2 = gsn_csm_xy(wks,x_axis(0:count_exp1(1)-1),track_error_avg_st_fh(1,i,0:count_exp1(1)-1),res) overlay(plot0,plot2) end if if (num_exp_cases.ge.3) then res@xyLineColor = colors_plot(3) res@xyLineThicknessF = 8.0 plot3 = gsn_csm_xy(wks,x_axis(0:count_exp1(2)-1),track_error_avg_st_fh(2,i,0:count_exp1(2)-1),res) overlay(plot0,plot3) end if if (num_exp_cases.ge.4) then res@xyLineColor = colors_plot(4) res@xyLineThicknessF = 8.0 plot4 = gsn_csm_xy(wks,x_axis(0:count_exp1(3)-1),track_error_avg_st_fh(3,i,0:count_exp1(3)-1),res) overlay(plot0,plot4) end if ; overlay 95% confidence interval polyres = True polyres@gsMarkerIndex = 1 polyres@gsMarkerSizeF = 0.04 bar_length = (/0.4,0.3,0.2,0.1/) ; Attach a legend marker_value = (/3,4,6,9/) mres = True txres = True txres@txFontHeightF = 0.016 do k=1,dimsizes(labels2)-1 ; start at one to skip CTL label txres@txFontColor = colors_plot(k) gsn_text_ndc(wks,labels2(k),xtxt3(k),ytxt3(k),txres) if (k.ge.1) then mres@gsMarkerIndex = marker_value(k-1) mres@gsMarkerSizeF = 16.0 mres@gsMarkerThicknessF = 4.0 mres@gsMarkerColor = colors_plot(k) gsn_polymarker_ndc(wks,xtxt3(k)-0.09,ytxt3(k),mres) end if end do do k=0,dimsizes(EXP_name)-1 ; loop over experiments polyres@gsLineColor = colors_plot(k+1) polyres@gsLineThicknessF = 4.0 do w=0,x_axis_max if (.not.ismissing(track_error_avg_st_fh(k,i,w))) then gsn_polyline(wks,plot0,(/x_axis(w),x_axis(w)/),(/1.0*interval_track_st_fh(k,i,w),-1.0*interval_track_st_fh(k,i,w)/),polyres) gsn_polyline(wks,plot0,(/x_axis(w)-bar_length(k),x_axis(w)+bar_length(k)/),(/1.0*interval_track_st_fh(k,i,w),1.0*interval_track_st_fh(k,i,w)/),polyres) ; error bar along x gsn_polyline(wks,plot0,(/x_axis(w)-bar_length(k),x_axis(w)+bar_length(k)/),(/-1.0*interval_track_st_fh(k,i,w),-1.0*interval_track_st_fh(k,i,w)/),polyres) ; error bar along x end if end do end do draw(plot0) frame(wks) delete(res) delete(polyres) ; intensity error max_diff = max(intensity_error_avg_st_fh(:,i,:)) ; set bounds min_diff = min(intensity_error_avg_st_fh(:,i,:)) range_diff = max_diff - min_diff wks = gsn_open_wks (output_type,"intensity_confidence_interval_fcst_hour-avg-" + g5nr_name(i)) res = True ; plot mods desired res@gsnXYBarChart = False res@tiMainString = "Intensity Error Difference (" + g5nr_name(i) + ")" res@trYMaxF = max_diff + 1.2*range_diff res@trYMinF = min_diff - 1.2*range_diff res@gsnFrame = False ; don't advance frame yet res@gsnDraw = False res@gsnFrame = False res@vpXF = 0.15 res@vpYF = 0.90 res@vpWidthF = 0.7 res@vpHeightF = 0.6 res@gsnScale = True res@trXMinF = 0 res@trXMaxF = 28 x_axis = ispan(0,28,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour res@tmXBValues = ispan(0,28,2) res@tmXBLabels = ispan(0,168,12) res@tiYAxisString = "Intensity Error Difference (hPa)" res@tiXAxisString = "Forecast lead time (hr)" res@tiXAxisFontHeightF = 0.015 res@tmYROn = False res@tmYROn = False res@tmXMajorGrid = True ; implement x grid res@tmXMajorGridThicknessF = 1.0 ; 2.0 is default res@tmXMajorGridLineDashPattern = 2 res@tmYRMinorOn = False res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 8.0 res@tiMainFontHeightF = 0.021 res@xyLineColor = colors_plot(1) res@gsnYRefLine = (/0/) ; first plot plot0 = gsn_csm_xy(wks,x_axis(0:count_exp1(0)-1),intensity_error_avg_st_fh(0,i,0:count_exp1(0)-1),res) if (num_exp_cases.ge.2) then res@xyLineColor = colors_plot(2) res@xyLineThicknessF = 8.0 plot2 = gsn_csm_xy(wks,x_axis(0:count_exp1(1)-1),intensity_error_avg_st_fh(1,i,0:count_exp1(1)-1),res) overlay(plot0,plot2) end if if (num_exp_cases.ge.3) then res@xyLineColor = colors_plot(3) res@xyLineThicknessF = 8.0 plot3 = gsn_csm_xy(wks,x_axis(0:count_exp1(2)-1),intensity_error_avg_st_fh(2,i,0:count_exp1(2)-1),res) overlay(plot0,plot3) end if if (num_exp_cases.ge.4) then res@xyLineColor = colors_plot(4) res@xyLineThicknessF = 8.0 plot4 = gsn_csm_xy(wks,x_axis(0:count_exp1(3)-1),intensity_error_avg_st_fh(3,i,0:count_exp1(3)-1),res) overlay(plot0,plot4) end if ; overlay 95% confidence interval polyres = True polyres@gsMarkerIndex = 1 polyres@gsMarkerSizeF = 0.04 bar_length = (/0.4,0.3,0.2,0.1/) ; Attach a legend marker_value = (/3,4,6,9/) mres = True txres = True txres@txFontHeightF = 0.016 do k=1,dimsizes(labels2)-1 ; start at one to skip CTL label txres@txFontColor = colors_plot(k) gsn_text_ndc(wks,labels2(k),xtxt3(k),ytxt3(k),txres) if (k.ge.1) then mres@gsMarkerIndex = marker_value(k-1) mres@gsMarkerSizeF = 16.0 mres@gsMarkerThicknessF = 4.0 mres@gsMarkerColor = colors_plot(k) gsn_polymarker_ndc(wks,xtxt3(k)-0.09,ytxt3(k),mres) end if end do do k=0,dimsizes(EXP_name)-1 ; loop over experiments polyres@gsLineColor = colors_plot(k+1) polyres@gsLineThicknessF = 4.0 do w=0,x_axis_max if (.not.ismissing(intensity_error_avg_st_fh(k,i,w))) then gsn_polyline(wks,plot0,(/x_axis(w),x_axis(w)/),(/1.0*interval_intensity_st_fh(k,i,w),-1.0*interval_intensity_st_fh(k,i,w)/),polyres) gsn_polyline(wks,plot0,(/x_axis(w)-bar_length(k),x_axis(w)+bar_length(k)/),(/1.0*interval_intensity_st_fh(k,i,w),1.0*interval_intensity_st_fh(k,i,w)/),polyres) ; error bar along x gsn_polyline(wks,plot0,(/x_axis(w)-bar_length(k),x_axis(w)+bar_length(k)/),(/-1.0*interval_intensity_st_fh(k,i,w),-1.0*interval_intensity_st_fh(k,i,w)/),polyres) ; error bar along x end if end do end do draw(plot0) frame(wks) delete(res) delete(polyres) ; wind error max_diff = max(wind_error_avg_st_fh(:,i,:)) ; set bounds min_diff = min(wind_error_avg_st_fh(:,i,:)) range_diff = max_diff - min_diff wks = gsn_open_wks (output_type,"wind_confidence_interval_fcst_hour-avg-" + g5nr_name(i)) res = True ; plot mods desired res@gsnXYBarChart = False res@tiMainString = "Wind Error Difference (" + g5nr_name(i) + ")" res@trYMaxF = max_diff + 1.2*range_diff res@trYMinF = min_diff - 1.2*range_diff res@gsnFrame = False ; don't advance frame yet res@gsnDraw = False res@gsnFrame = False res@vpXF = 0.15 res@vpYF = 0.90 res@vpWidthF = 0.7 res@vpHeightF = 0.6 res@gsnScale = True res@trXMinF = 0 res@trXMaxF = 28 x_axis = ispan(0,28,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour res@tmXBValues = ispan(0,28,2) res@tmXBLabels = ispan(0,168,12) res@tiYAxisString = "Wind Error Difference (kts)" res@tiXAxisString = "Forecast lead time (hr)" res@tiXAxisFontHeightF = 0.015 res@tmYROn = False res@tmYROn = False res@tmXMajorGrid = True ; implement x grid res@tmXMajorGridThicknessF = 1.0 ; 2.0 is default res@tmXMajorGridLineDashPattern = 2 res@tmYRMinorOn = False res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 8.0 res@tiMainFontHeightF = 0.021 res@xyLineColor = colors_plot(1) res@gsnYRefLine = (/0/) ; first plot plot0 = gsn_csm_xy(wks,x_axis(0:count_exp1(0)-1),wind_error_avg_st_fh(0,i,0:count_exp1(0)-1),res) if (num_exp_cases.ge.2) then res@xyLineColor = colors_plot(2) res@xyLineThicknessF = 8.0 plot2 = gsn_csm_xy(wks,x_axis(0:count_exp1(1)-1),wind_error_avg_st_fh(1,i,0:count_exp1(1)-1),res) overlay(plot0,plot2) end if if (num_exp_cases.ge.3) then res@xyLineColor = colors_plot(3) res@xyLineThicknessF = 8.0 plot3 = gsn_csm_xy(wks,x_axis(0:count_exp1(2)-1),wind_error_avg_st_fh(2,i,0:count_exp1(2)-1),res) overlay(plot0,plot3) end if if (num_exp_cases.ge.4) then res@xyLineColor = colors_plot(4) res@xyLineThicknessF = 8.0 plot4 = gsn_csm_xy(wks,x_axis(0:count_exp1(3)-1),wind_error_avg_st_fh(3,i,0:count_exp1(3)-1),res) overlay(plot0,plot4) end if ; Attach a legend marker_value = (/3,4,6,9/) mres = True txres = True txres@txFontHeightF = 0.016 do k=1,dimsizes(labels2)-1 ; start at one to skip CTL label txres@txFontColor = colors_plot(k) gsn_text_ndc(wks,labels2(k),xtxt3(k),ytxt3(k),txres) if (k.ge.1) then mres@gsMarkerIndex = marker_value(k-1) mres@gsMarkerSizeF = 16.0 mres@gsMarkerThicknessF = 4.0 mres@gsMarkerColor = colors_plot(k) gsn_polymarker_ndc(wks,xtxt3(k)-0.09,ytxt3(k),mres) end if end do ; overlay 95% confidence interval polyres = True polyres@gsMarkerIndex = 1 polyres@gsMarkerSizeF = 0.04 bar_length = (/0.4,0.3,0.2,0.1/) do k=0,dimsizes(EXP_name)-1 ; loop over experiments polyres@gsLineColor = colors_plot(k+1) polyres@gsLineThicknessF = 4.0 do w=0,x_axis_max if (.not.ismissing(wind_error_avg_st_fh(k,i,w))) then gsn_polyline(wks,plot0,(/x_axis(w),x_axis(w)/),(/1.0*interval_wind_st_fh(k,i,w),-1.0*interval_wind_st_fh(k,i,w)/),polyres) gsn_polyline(wks,plot0,(/x_axis(w)-bar_length(k),x_axis(w)+bar_length(k)/),(/1.0*interval_wind_st_fh(k,i,w),1.0*interval_wind_st_fh(k,i,w)/),polyres) ; error bar along x gsn_polyline(wks,plot0,(/x_axis(w)-bar_length(k),x_axis(w)+bar_length(k)/),(/-1.0*interval_wind_st_fh(k,i,w),-1.0*interval_wind_st_fh(k,i,w)/),polyres) ; error bar along x end if end do end do draw(plot0) frame(wks) delete(res) delete(polyres) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; delete([/track_temp,max_track,min_track,max_wind,min_wind,wind_temp,max_int,min_int,int_temp,max_fcst,min_fcst,x_axis/]) end if end if delete(location_max_st_temp) end do print("-------------------------------------------------------") print("Total storms analyzed: " + count_storms) print("-------Track Program completed-------------------------") end