; Reads the best storm tracks of tropical cyclones from G5NR ; and from the GFS experiment runs ; and computes an average track, wind, and intensity error ; as well as individual results for each case ; ; PRIOR TO RUNNING THIS CODE, RUN THIS COMMAND ON YOUR ATCFUNIX.GFS.* FILES TO GET ALL THE TRACKS INTO ONE TEXT FILE ; "grep "AL, 01" atcfunix.gfs.* > AL01.txt, AL01 refers to the storm name in your runs, will need to know this... ; ; NOTE: THIS PROGRAM DOES NOT ACCOUNT FOR THE FACT THAT THE TRACK FILE (TCVITALS PRODUCED ABOVE) IN GFS HAS DUPLICATE FORECAST LINES ; TO FIX THIS: USE LINUX COMMAND: sort -u -t',' -k3,6 WP13.txt > WP13_fix.txt AS AN EXAMPLE TO CREATE NEW ; FILE. DO THIS BEFORE RUNNING PROGRAM AND AFTER YOU DO THE GREP COMMAND ; ; IF YOU WANT TO RUN PROGRAM FOR SPECIFIC BASINS, PLACE TEXT FILES IN THEIR OWN DIRECTORY AS INDICATED BELOW ; AND RUN ACCORDINGLY ; ; RESULTS ARE PRESENTED OUT TO 168 HOURS SINCE THE LONGEST LENGTH IN THE TCVITALS IS 168 HRS ; ; IF FUTURE TCVITALS FILES OR G5NR BEST TRACK FILE FORMAT IS CHANGED, PROGRAM WILL NEED TO BE ADJUSTED ; ; ; Andrew Kren ; September 13, 2017 ; This version of the program can work on more than 1 experiment case at once, but still only one control begin ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; begin user's settings ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; if running by basin, change paths here to reflect which paths to run on ; make sure to only have the appropriate text files in the basin directory region_name = "All Basins" ; not case sensitive but make it easy ; set as "All Basins", "Eastern Pacific", "Atlantic", or "Western Pacific". Used for the plot titles, legends, etc. region_path = "ALL" ; CASE SENSITIVE ; for Atlantic -------- use ATL ; for Eastern Pacific - use EPA ; for Western Paficic - use WEP ; for all basins ------ use ALL CTL_name = "CTL" ; name for your control case, used for plot labels, etc., suggest to keep same name as folder below for consistency ; modify paths accordingly for location of track files path_g5nr = "/scratch4/BMC/qosap/Andrew.Kren/ncl_scripts/G5NR/" + region_path + "/" path_ctl = "/scratch4/BMC/qosap/Andrew.Kren/ncl_scripts/CTL/" + region_path + "/" ; experiment settings num_exp_cases = 2 ; number of exper not counting the control experiment EXP_name = (/"RO_PER","RO_ERR"/) ; name and path name for your experiment(s), if adding more than one, use quotes to separate each one path_exper = "/scratch4/BMC/qosap/Andrew.Kren/ncl_scripts/" ; location of output files, make sure all the folders for appropriate cases are inside this subdirectory output_type = "png" ; "x11", "ps", "eps", "pdf", etc. ; format of the output images you want ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; end of user's settings ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; begin of main program ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; list files in directory to determine how many dropsondes to read fils_g5nr = systemfunc("csh -c 'cd " + path_g5nr + " ; ls *.dat'") sizes_g5nr = dimsizes(fils_g5nr) ; number of g5nr storms determined from the files above fils_ctl = systemfunc("csh -c 'cd " + path_ctl + " ; ls *fix*.txt'") ; extenstion *fix* since those are altered based on linux command sort above sizes_ctl = dimsizes(fils_ctl) region_file = str_join( str_split(region_name," "),"_") ; output file name for images ; determine how many experiment cases not counting ctl fils_exp = new((/dimsizes(EXP_name),sizes_ctl/),"string") path_exp = new((/dimsizes(EXP_name)/),"string") input_file_exp = fils_exp do a=0,num_exp_cases-1 ; loop over experiment cases to determine all fils_exp path_exp(a) = path_exper + EXP_name(a) + "/" + region_path + "/" temp = systemfunc("csh -c 'cd " + path_exp(a) + " ; ls *fix*.txt'") ; original files are the non-fix named files in the same directory fils_exp(a,0:dimsizes(temp)-1) = temp input_file_exp(a,0:dimsizes(temp)-1) = path_exp(a) + fils_exp(a,0:dimsizes(temp)-1) delete(temp) end do xtxt3 = (/0.35,0.35,0.35,0.35,0.35,0.35,0.35/) ytxt3 = (/0.77,0.74,0.71,0.68,0.65,0.62,0.59/) input_file_g5nr = path_g5nr + fils_g5nr input_file_ctl = path_ctl + fils_ctl scale_factor = 0.1 scale_factor_west = -0.1 fcst_hr = ispan(0,168,6) ; indices 0 to 28 ; forecast hours for program ; create arrays to store storm track data g5nr_pres = new((/sizes_g5nr,300/),"float") g5nr_name = new((/sizes_g5nr/),"string") g5nr_lat = g5nr_pres g5nr_lon = g5nr_pres g5nr_wind = g5nr_pres g5nr_ymd = new((/sizes_g5nr,300/),"string") ctl_pres = new((/sizes_g5nr,3000/),"float") ctl_name = new((/sizes_g5nr/),"string") ctl_lat = ctl_pres ctl_lon = ctl_pres ctl_wind = ctl_pres ctl_ymd = new((/sizes_g5nr,3000/),"string") ctl_hr = ctl_ymd exp_pres = new((/dimsizes(EXP_name),sizes_g5nr,3000/),"float") exp_name = new((/dimsizes(EXP_name),sizes_g5nr/),"string") exp_lat = exp_pres exp_lon = exp_pres exp_wind = exp_pres exp_ymd = new((/dimsizes(EXP_name),sizes_g5nr,3000/),"string") exp_hr = exp_ymd track_error_ctl = new((/sizes_g5nr,20,29/),"float") ; saved track errors and intensity errors, 29 corresponds to 0 to 168 hours of 6hr output intensity_error_ctl = track_error_ctl ; array structure of (storm,cycle,forcast_hr) track_error_exp = new((/dimsizes(EXP_name),sizes_g5nr,20,29/),"float") ctl_ymd_array = new((/sizes_g5nr,20/),"string") exp_ymd_array = new((/dimsizes(EXP_name),sizes_g5nr,20/),"string") fcst_hr_array_ctl = new((/29/),"float") fcst_hr_array_exp = new((/dimsizes(EXP_name),29/),"float") fcst_hr_array_ctl2 = fcst_hr_array_ctl fcst_hr_array_exp2 = fcst_hr_array_exp intensity_error_exp = track_error_exp wind_error_exp = track_error_exp wind_error_ctl = track_error_ctl ctl_lat_array = track_error_ctl ctl_lon_array = track_error_ctl exp_lat_array = track_error_exp exp_lon_array = track_error_exp ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; loop over each of the storms, plotting the storm track at each iteration ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; print("Reading G5NR storm tracks...") do i=0,dimsizes(fils_g5nr)-1 ; read in data for this storm data1 = asciiread(input_file_g5nr(i),-1,"string") ; comma separated values delim = "," g5nr_ymd(i,0:dimsizes(data1)-1) = str_get_field(data1,3,delim) g5nr_pres(i,0:dimsizes(data1)-1) = tofloat(str_get_field(data1,10,delim)) g5nr_wind(i,0:dimsizes(data1)-1) = tofloat(str_get_field(data1,9,delim)) ; knots name_1 = str_squeeze(str_get_field(data1,1,delim)) name_2 = str_squeeze(str_get_field(data1,2,delim)) g5nr_name(i) = str_squeeze(name_1(0)+name_2(0)) g5nr_lat(i,0:dimsizes(data1)-1) = tofloat(str_get_cols(str_get_field(data1,7,delim),0,3)) * scale_factor lon = str_get_cols(str_get_field(data1,8,delim),0,4) deg_we = str_get_cols(str_get_field(data1,8,delim),5,5) do j=0,dimsizes(data1)-1 ; assign correct scale factor to longitudes based on degW or degE if (deg_we(j).eq."E") then lon(j) = tofloat(lon(j)) * scale_factor ; deg W end if if (deg_we(j).eq."W") then lon(j) = tofloat(lon(j)) * scale_factor_west ; deg E end if end do g5nr_lon(i,0:dimsizes(data1)-1) = tofloat(lon) delete([/data1,lon,deg_we/]) delete([/name_1,name_2/]) end do ; end loop over G5NR storms ; read the control and experiment storm track files ; CTL case reading print("Reading CTL storm tracks...") do i=0,dimsizes(fils_ctl)-1 ; open file and read data1 = asciiread(input_file_ctl(i),-1,"string") ; 20 column array of values delim = "," ctl_ymd(i,0:dimsizes(data1)-1) = str_get_field(data1,3,delim) ctl_hr(i,0:dimsizes(data1)-1) = str_get_field(data1,6,delim) ctl_lat(i,0:dimsizes(data1)-1) = tofloat(str_get_cols(str_get_field(data1,7,delim),0,3)) * scale_factor lon = str_get_cols(str_get_field(data1,8,delim),0,4) ctl_pres(i,0:dimsizes(data1)-1) = tofloat(str_get_field(data1,10,delim)) ctl_wind(i,0:dimsizes(data1)-1) = tofloat(str_get_field(data1,9,delim)) ; knots name_1 = str_squeeze(str_get_cols(str_get_field(data1,1,delim),24,25)) name_2 = str_squeeze(str_get_field(data1,2,delim)) ctl_name(i) = str_squeeze(name_1(0)+name_2(0)) delete([/name_1,name_2/]) deg_we = str_get_cols(str_get_field(data1,8,delim),5,5) do j=0,dimsizes(data1)-1 ; assign correct scale factor to longitudes based on degW or degE if (deg_we(j).eq."E") then lon(j) = tofloat(lon(j)) * scale_factor ; deg W end if if (deg_we(j).eq."W") then lon(j) = tofloat(lon(j)) * scale_factor_west ; deg E end if end do ctl_lon(i,0:dimsizes(data1)-1) = tofloat(lon) delete([/data1,lon,deg_we/]) end do ; EXP cases reading print("Reading EXP(s) storm tracks...") do a=0,dimsizes(EXP_name)-1 do i=0,dimsizes(fils_exp(a,:))-1 ; loop over storms for this experiment case ; open file and read data1 = asciiread(input_file_exp(a,i),-1,"string") ; 20 column array of values delim = "," exp_ymd(a,i,0:dimsizes(data1)-1) = str_get_field(data1,3,delim) exp_hr(a,i,0:dimsizes(data1)-1) = str_get_field(data1,6,delim) exp_lat(a,i,0:dimsizes(data1)-1) = tofloat(str_get_cols(str_get_field(data1,7,delim),0,3)) * scale_factor lon = str_get_cols(str_get_field(data1,8,delim),0,4) exp_pres(a,i,0:dimsizes(data1)-1) = tofloat(str_get_field(data1,10,delim)) exp_wind(a,i,0:dimsizes(data1)-1) = tofloat(str_get_field(data1,9,delim)) ; knots name_1 = str_squeeze(str_get_cols(str_get_field(data1,1,delim),24,25)) name_2 = str_squeeze(str_get_field(data1,2,delim)) exp_name(a,i) = str_squeeze(name_1(0)+name_2(0)) delete([/name_1,name_2/]) deg_we = str_get_cols(str_get_field(data1,8,delim),5,5) do j=0,dimsizes(data1)-1 ; assign correct scale factor to longitudes based on degW or degE if (deg_we(j).eq."E") then lon(j) = tofloat(lon(j)) * scale_factor ; deg W end if if (deg_we(j).eq."W") then lon(j) = tofloat(lon(j)) * scale_factor_west ; deg E end if end do exp_lon(a,i,0:dimsizes(data1)-1) = tofloat(lon) delete([/data1,lon,deg_we/]) end do ; end loop over each experiment end do ; end loop over experiment forecasts ;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; 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 ;;;;;;;;;;;;;;;;;;;;;;;; ;; 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)/),"float") 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 = min(count_exp1) check_check = check_exp 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@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 plot1 = gsn_csm_xy(wks,x_axis,track_error_exp(0,i,check_exp(0),0:x_axis_max),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(1,i,check_exp(1),0:x_axis_max),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(2,i,check_exp(2),0:x_axis_max),res) overlay(plot0,plot3) end if if (num_exp_cases.ge.4) then res@xyLineColor = colors_plot(4) res@xyLineThicknessF = 8.0 plot2 = gsn_csm_xy(wks,x_axis,track_error_exp(3,i,check_exp(3),0:x_axis_max),res) 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) 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 plot1 = gsn_csm_xy(wks,x_axis,wind_error_exp(0,i,check_exp(0),0:x_axis_max),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(1,i,check_exp(1),0:x_axis_max),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(2,i,check_exp(2),0:x_axis_max),res) overlay(plot0,plot3) end if if (num_exp_cases.ge.4) then res@xyLineColor = colors_plot(4) res@xyLineThicknessF = 8.0 plot2 = gsn_csm_xy(wks,x_axis,wind_error_exp(3,i,check_exp(3),0:x_axis_max),res) 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@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 plot1 = gsn_csm_xy(wks,x_axis,intensity_error_exp(0,i,check_exp(0),0:x_axis_max),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(1,i,check_exp(1),0:x_axis_max),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(2,i,check_exp(2),0:x_axis_max),res) overlay(plot0,plot3) end if if (num_exp_cases.ge.4) then res@xyLineColor = colors_plot(4) res@xyLineThicknessF = 8.0 plot2 = gsn_csm_xy(wks,x_axis,intensity_error_exp(3,i,check_exp(3),0:x_axis_max),res) 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:length_plot_ctl),exp_lat_array(0,i,check_exp(0),0:length_plot_ctl),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:length_plot_ctl),exp_lat_array(1,i,check_exp(1),0:length_plot_ctl),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:length_plot_ctl),exp_lat_array(2,i,check_exp(2),0:length_plot_ctl),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:length_plot_ctl),exp_lat_array(3,i,check_exp(3),0:length_plot_ctl),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:length_plot_ctl),exp_lat_array(0,i,check_exp(0),0:length_plot_ctl),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:length_plot_ctl),exp_lat_array(1,i,check_exp(1),0:length_plot_ctl),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:length_plot_ctl),exp_lat_array(2,i,check_exp(2),0:length_plot_ctl),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:length_plot_ctl),exp_lat_array(3,i,check_exp(3),0:length_plot_ctl),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...") track_error_exp!0 = "exp" track_error_exp!1 = "storm" track_error_exp!2 = "cycle" track_error_exp!3 = "fcst_hr" intensity_error_exp!0 = "exp" intensity_error_exp!1 = "storm" intensity_error_exp!2 = "cycle" intensity_error_exp!3 = "fcst_hr" wind_error_exp!0 = "exp" wind_error_exp!1 = "storm" wind_error_exp!2 = "cycle" wind_error_exp!3 = "fcst_hr" track_error_exp_reorder = track_error_exp(storm|:,fcst_hr|:,exp|:,cycle|:) intensity_error_exp_reorder = intensity_error_exp(storm|:,fcst_hr|:,exp|:,cycle|:) wind_error_exp_reorder = wind_error_exp(storm|:,fcst_hr|:,exp|:,cycle|:) track_error_ctl_st_fh = dim_avg_n(track_error_ctl,1) ; 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,1) intensity_error_exp_st_fh = dim_avg_n(intensity_error_exp_reorder,(/3/)) wind_error_ctl_st_fh = dim_avg_n(wind_error_ctl,1) wind_error_exp_st_fh = dim_avg_n(wind_error_exp_reorder,(/3/)) track_error_ctl_fh = dim_avg_n(track_error_ctl,(/0,1/)) ; average as function of forecast hour only track_error_exp_fh = dim_avg_n(track_error_exp,(/1,2/)) track_error_ctl_fh_var = dim_variance_n(track_error_ctl,(/0,1/)) ; average as function of forecast hour only track_error_exp_fh_var = dim_variance_n(track_error_exp,(/1,2/)) track_error_ctl_fh_std = dim_stddev_n(track_error_ctl,(/0,1/)) ; average as function of forecast hour only track_error_exp_fh_std = dim_stddev_n(track_error_exp,(/1,2/)) intensity_error_ctl_fh_var = dim_variance_n(intensity_error_ctl,(/0,1/)) ; average as function of forecast hour only intensity_error_exp_fh_var = dim_variance_n(intensity_error_exp,(/1,2/)) intensity_error_ctl_fh_std = dim_stddev_n(intensity_error_ctl,(/0,1/)) ; average as function of forecast hour only intensity_error_exp_fh_std = dim_stddev_n(intensity_error_exp,(/1,2/)) intensity_error_ctl_fh = dim_avg_n(intensity_error_ctl,(/0,1/)) intensity_error_exp_fh = dim_avg_n(intensity_error_exp,(/1,2/)) wind_error_ctl_fh_var = dim_variance_n(wind_error_ctl,(/0,1/)) ; average as function of forecast hour only wind_error_exp_fh_var = dim_variance_n(wind_error_exp,(/1,2/)) wind_error_ctl_fh_std = dim_stddev_n(wind_error_ctl,(/0,1/)) ; average as function of forecast hour only wind_error_exp_fh_std = dim_stddev_n(wind_error_exp,(/1,2/)) wind_error_ctl_fh = dim_avg_n(wind_error_ctl,(/0,1/)) wind_error_exp_fh = dim_avg_n(wind_error_exp,(/1,2/)) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; total average over all storms and hours, one value here track_error_ctl_avg = avg(track_error_ctl) ; total average track_error_exp_avg = dim_avg_n(track_error_exp,(/1,2,3/)) intensity_error_ctl_avg = avg(intensity_error_ctl) intensity_error_exp_avg = dim_avg_n(intensity_error_exp,(/1,2,3/)) wind_error_ctl_avg = avg(wind_error_ctl) wind_error_exp_avg = dim_avg_n(wind_error_exp,(/1,2,3/)) ; determine number of non-missing values at each forecast hour in each experiment do i=0,28 ; loop over forecast hours (output every 6 hours) fcst_hr_array_ctl(i) = num(.not.ismissing(track_error_ctl(:,:,i))) end do 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(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 siglvl = 0.05 ; 95% significance level do a=0,dimsizes(EXP_name)-1 track_error_diff_std(a,:) = dim_stddev_n(track_error_exp(a,:,:,:)-track_error_ctl,(/0,1/)) ; forecast hour only intensity_error_diff_std(a,:) = dim_stddev_n(intensity_error_exp(a,:,:,:)-intensity_error_ctl,(/0,1/)) ; forecast hour only, wind_error_diff_std(a,:) = dim_stddev_n(wind_error_exp(a,:,:,:)-wind_error_ctl,(/0,1/)) ; forecast hour only track_error_diff_avg(a,:) = dim_avg_n(track_error_exp(a,:,:,:)-track_error_ctl,(/0,1/)) ; forecast hour only intensity_error_diff_avg(a,:) = dim_avg_n(intensity_error_exp(a,:,:,:)-intensity_error_ctl,(/0,1/)) ; forecast hour only wind_error_diff_avg(a,:) = dim_avg_n(wind_error_exp(a,:,:,:)-wind_error_ctl,(/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(y).le.fcst_hr_array_exp(a,y)) then number_sample(a,y) = fcst_hr_array_ctl(y) else number_sample(a,y) = fcst_hr_array_exp(a,y) end if end do interval_track_diff(a,:) = 2.0*track_error_diff_std(a,:) / sqrt(number_sample(a,:)) ; % 95 interval interval_intensity_diff(a,:) = 2.0*intensity_error_diff_std(a,:) / sqrt(number_sample(a,:)) ; 95% interval interval_wind_diff(a,:) = 2.0*wind_error_diff_std(a,:) / sqrt(number_sample(a,:)) ; 95% interval print(track_error_diff_avg(a,:)) track_error_diff_avg(a,:) = track_error_exp_fh(a,:)-track_error_ctl_fh print(track_error_diff_avg(a,:)) exit intensity_error_diff_avg(a,:) = abs(intensity_error_exp_fh(a,:)-intensity_error_ctl_fh) wind_error_diff_avg(a,:) = abs(wind_error_exp_fh(a,:)-wind_error_ctl_fh) ave_ctl = track_error_ctl_fh(0:28) ave_exp = track_error_exp_fh(a,0:28) var_ctl = track_error_ctl_fh_var(0:28) var_exp = track_error_exp_fh_var(a,0:28) s_ctl = fcst_hr_array_ctl(0:28) s_exp = fcst_hr_array_exp(a,0:28) iflag = False tval_opt = False prob_fh_track(a,:) = ttest(ave_ctl,var_ctl,s_ctl,ave_exp,var_exp,s_exp,iflag,tval_opt) ;print("Significance p-value over all storms for track over forecast hour is: " + prob_fh_track) delete([/s_ctl,s_exp,ave_ctl,ave_exp,var_ctl,var_exp/]) ; intensity ; compute significance over full avg and over forecast hour for track and intensity ave_ctl = intensity_error_ctl_fh(0:28) ave_exp = intensity_error_exp_fh(a,0:28) var_ctl = intensity_error_ctl_fh_var(0:28) var_exp = intensity_error_exp_fh_var(a,0:28) s_ctl = fcst_hr_array_ctl(0:28) s_exp = fcst_hr_array_exp(a,0:28) iflag = False tval_opt = False prob_fh_intensity(a,:) = ttest(ave_ctl,var_ctl,s_ctl,ave_exp,var_exp,s_exp,iflag,tval_opt) ;print("Significance p-value over all storms for intensity over forecast hour is: " + prob_fh_intensity) delete([/s_ctl,s_exp,ave_ctl,ave_exp,var_ctl,var_exp/]) ; wind ; compute significance over full avg and over forecast hour for track and intensity ave_ctl = wind_error_ctl_fh(0:28) ave_exp = wind_error_exp_fh(a,0:28) var_ctl = wind_error_ctl_fh_var(0:28) var_exp = wind_error_exp_fh_var(a,0:28) s_ctl = fcst_hr_array_ctl(0:28) s_exp = fcst_hr_array_exp(a,0:28) iflag = False tval_opt = False prob_fh_wind(a,:) = ttest(ave_ctl,var_ctl,s_ctl,ave_exp,var_exp,s_exp,iflag,tval_opt) ;print("Significance p-value over all storms for intensity over forecast hour is: " + prob_fh_intensity) delete([/s_ctl,s_exp,ave_ctl,ave_exp,var_ctl,var_exp/]) 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,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,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,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@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@tiMainFontHeightF = 0.021 res@xyLineColor = colors_plot(0) plot0 = gsn_csm_xy(wks,x_axis,track_error_ctl_fh(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@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@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 = "# Samples" bres@tiYAxisSide = "Right" bres@tiYAxisAngleF = 270 bres@tmXTOn = False plot5 = gsn_csm_xy(wks,x_axis,fcst_hr_array_ctl(0:28),bres) bres@gsnXYBarChartColors = "firebrick1" ;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 = "orange" ; 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 = "limegreen" ; 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 = "lightblue" ; 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 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) end do ; attach significance marker_value = (/3,4,6,9/) do k=0,dimsizes(EXP_name)-1 do w=0,28 ; loop over forecast hours mres = True mres@gsMarkerIndex = marker_value(k) mres@gsMarkerSizeF = 16.0 mres@gsMarkerThicknessF = 4.0 mres@gsMarkerColor = colors_plot(k+1) ;if (.not.ismissing(prob_fh_track(k,w)) .and. prob_fh_track(k,w).le.siglvl) then ; then label as significant ; gsn_polymarker(wks,plot0,x_axis(w),20.0,mres) ;end if if (.not.ismissing(track_error_diff_avg(k,w)) .and. track_error_diff_avg(k,w).gt.interval_track_diff(k,w)) then ; 95% confidence interval gsn_polymarker(wks,plot0,x_axis(w),20.0,mres) end if end do end do draw(plot0) frame(wks) delete(res) delete(bres) ; 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@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@xyLineThicknessF = 8.0 res@tiMainFontHeightF = 0.021 res@xyLineColor = colors_plot(0) plot0 = gsn_csm_xy(wks,x_axis,wind_error_ctl_fh(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@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@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 = "# Samples" bres@tiYAxisSide = "Right" bres@tiYAxisAngleF = 270 bres@tmXTOn = False plot5 = gsn_csm_xy(wks,x_axis,fcst_hr_array_ctl(0:28),bres) bres@gsnXYBarChartColors = "firebrick1" 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 = "orange" 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 = "limegreen" 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 = "lightblue" 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 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) end do ; attach significance marker_value = (/3,4,6,9/) do k=0,dimsizes(EXP_name)-1 do w=0,28 ; loop over forecast hours mres = True mres@gsMarkerIndex = marker_value(k) mres@gsMarkerSizeF = 16.0 mres@gsMarkerThicknessF = 4.0 mres@gsMarkerColor = colors_plot(k+1); blue3 ;if (.not.ismissing(prob_fh_wind(k,w)) .and. prob_fh_wind(k,w).le.siglvl) then ; then label as significant ; gsn_polymarker(wks,plot0,x_axis(w),20.0,mres) ;end if if (.not.ismissing(wind_error_diff_avg(k,w)) .and. wind_error_diff_avg(k,w).gt.interval_wind_diff(k,w)) then ; 95% confidence interval gsn_polymarker(wks,plot0,x_axis(w),5.0,mres) end if end do end do draw(plot0) frame(wks) delete(res) delete(bres) ; 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 = "Track 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 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@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(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@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@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 = "# Samples" bres@tiYAxisSide = "Right" bres@tiYAxisAngleF = 270 bres@tmXTOn = False plot5 = gsn_csm_xy(wks,x_axis,fcst_hr_array_ctl(0:28),bres) bres@gsnXYBarChartColors = "firebrick1" 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 = "orange" 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 = "limegreen" 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 = "lightblue" 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 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) end do ; attach significance marker_value = (/3,4,6,9/) do k=0,dimsizes(EXP_name)-1 do w=0,28 ; loop over forecast hours mres = True mres@gsMarkerIndex = marker_value(k) mres@gsMarkerSizeF = 16.0 mres@gsMarkerThicknessF = 4.0 mres@gsMarkerColor = colors_plot(k+1) ;if (.not.ismissing(prob_fh_intensity(k,w)) .and. prob_fh_intensity(k,w).le.siglvl) then ; then label as significant ; gsn_polymarker(wks,plot0,x_axis(w),20.0,mres) ;end if if (.not.ismissing(intensity_error_diff_avg(k,w)) .and. intensity_error_diff_avg(k,w).gt.interval_intensity_diff(k,w)) then ; 95% confidence interval gsn_polymarker(wks,plot0,x_axis(w),5.0,mres) end if end do end do draw(plot0) frame(wks) delete(res) delete(bres) ; 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@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@tiMainFontHeightF = 0.021 res@xyLineColor = colors_plot(0) plot0 = gsn_csm_xy(wks,x_axis,track_error_ctl_fh(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)) max_std = (track_error_ctl_fh_std(ind_valid)) / sqrt(fcst_hr_array_ctl(ind_valid)) spread = (/track_error_ctl_fh(ind_valid)-max_std,track_error_ctl_fh(ind_valid)+max_std/) res@gsnXYFillOpacities = 0.45 ctl_spread = gsn_csm_xy(wks,x_axis(ind_valid),spread,res) res@gsnXYFillColors = "firebrick1" 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 = "orange" 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 = "limegreen" 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 = "lightblue" 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 = "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@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 = "# Samples" bres@tiYAxisSide = "Right" bres@tiYAxisAngleF = 270 bres@tmXTOn = False ;bres@gsnXYBarChartFillOpacityF = 0.2 plot5 = gsn_csm_xy(wks,x_axis,fcst_hr_array_ctl(0:28),bres) bres@gsnXYBarChartColors = "firebrick1" 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 = "orange" 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 = "limegreen" 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 = "lightblue" plot9 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp(3,0:28),bres) draw(plot9) 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),xtxt3(k),ytxt3(k),txres) end do ; attach significance marker_value = (/3,4,6,9/) do k=0,dimsizes(EXP_name)-1 do w=0,28 ; loop over forecast hours mres = True mres@gsMarkerIndex = marker_value(k) mres@gsMarkerSizeF = 16.0 mres@gsMarkerThicknessF = 4.0 mres@gsMarkerColor = colors_plot(k+1) ;if (.not.ismissing(prob_fh_track(k,w)) .and. prob_fh_track(k,w).le.siglvl) then ; then label as significant ; gsn_polymarker(wks,plot0,x_axis(w),20.0,mres) ;end if if (.not.ismissing(track_error_diff_avg(k,w)) .and. track_error_diff_avg(k,w).gt.interval_track_diff(k,w)) then ; 95% confidence interval gsn_polymarker(wks,plot0,x_axis(w),20.0,mres) end if end do end do draw(plot0) frame(wks) delete(res) delete(bres) ; 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@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@xyLineThicknessF = 8.0 res@tiMainFontHeightF = 0.021 res@xyLineColor = colors_plot(0) plot0 = gsn_csm_xy(wks,x_axis,wind_error_ctl_fh(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)) max_std = (wind_error_ctl_fh_std(ind_valid)) / sqrt(fcst_hr_array_ctl(ind_valid)) spread = (/wind_error_ctl_fh(ind_valid)-max_std,wind_error_ctl_fh(ind_valid)+max_std/) res@gsnXYFillOpacities = 0.45 ctl_spread = gsn_csm_xy(wks,x_axis(ind_valid),spread,res) res@gsnXYFillColors = "firebrick1" 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 = "orange" 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 = "limegreen" 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 = "lightblue" 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@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@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 = "# Samples" bres@tiYAxisSide = "Right" bres@tiYAxisAngleF = 270 bres@tmXTOn = False plot5 = gsn_csm_xy(wks,x_axis,fcst_hr_array_ctl(0:28),bres) bres@gsnXYBarChartColors = "firebrick1" 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 = "orange" 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 = "limegreen" 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 = "lightblue" plot9 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp(3,0:28),bres) draw(plot9) 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),xtxt3(k),ytxt3(k),txres) end do ; attach significance marker_value = (/3,4,6,9/) do k=0,dimsizes(EXP_name)-1 do w=0,28 ; loop over forecast hours mres = True mres@gsMarkerIndex = marker_value(k) mres@gsMarkerSizeF = 16.0 mres@gsMarkerThicknessF = 4.0 mres@gsMarkerColor = colors_plot(k+1) ;if (.not.ismissing(prob_fh_wind(k,w)) .and. prob_fh_wind(k,w).le.siglvl) then ; then label as significant ; gsn_polymarker(wks,plot0,x_axis(w),20.0,mres) ;end if if (.not.ismissing(wind_error_diff_avg(k,w)) .and. wind_error_diff_avg(k,w).gt.interval_wind_diff(k,w)) then ; 95% confidence interval gsn_polymarker(wks,plot0,x_axis(w),5.0,mres) end if end do end do draw(plot0) frame(wks) delete(res) delete(bres) ; 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@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@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(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)) max_std = (intensity_error_ctl_fh_std(ind_valid)) / sqrt(fcst_hr_array_ctl(ind_valid)) spread = (/intensity_error_ctl_fh(ind_valid)-max_std,intensity_error_ctl_fh(ind_valid)+max_std/) res@gsnXYFillOpacities = 0.45 ctl_spread = gsn_csm_xy(wks,x_axis(ind_valid),spread,res) res@gsnXYFillColors = "firebrick1" 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 = "orange" 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 = "limegreen" 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 = "lightblue" 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@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@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 = "# Samples" bres@tiYAxisSide = "Right" bres@tiYAxisAngleF = 270 bres@tmXTOn = False plot5 = gsn_csm_xy(wks,x_axis,fcst_hr_array_ctl(0:28),bres) bres@gsnXYBarChartColors = "firebrick1" 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 = "orange" 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 = "limegreen" 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 = "lightblue" plot9 = gsn_csm_xy(wks,x_axis,fcst_hr_array_exp(3,0:28),bres) draw(plot9) 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),xtxt3(k),ytxt3(k),txres) end do ; attach significance marker_value = (/3,4,6,9/) do k=0,dimsizes(EXP_name)-1 do w=0,28 ; loop over forecast hours mres = True mres@gsMarkerIndex = marker_value(k) mres@gsMarkerSizeF = 16.0 mres@gsMarkerThicknessF = 4.0 mres@gsMarkerColor = colors_plot(k+1) ;if (.not.ismissing(prob_fh_intensity(k,w)) .and. prob_fh_intensity(k,w).le.siglvl) then ; then label as significant ; gsn_polymarker(wks,plot0,x_axis(w),20.0,mres) ;end if if (.not.ismissing(intensity_error_diff_avg(k,w)) .and. intensity_error_diff_avg(k,w).gt.interval_intensity_diff(k,w)) then ; 95% confidence interval gsn_polymarker(wks,plot0,x_axis(w),5.0,mres) end if end do end do draw(plot0) frame(wks) delete(res) delete(bres) 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 = check_exp ;;;;;;;;;;;;;;; ;;;; 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 ; 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(a,:,:,:)-track_error_ctl,1) ; forecast hour only intensity_error_std_st_fh(a,:,:) = dim_stddev_n(intensity_error_exp(a,:,:,:)-intensity_error_ctl,1) ; forecast hour only, wind_error_std_st_fh(a,:,:) = dim_stddev_n(wind_error_exp(a,:,:,:)-wind_error_ctl,1) ; forecast hour only track_error_avg_st_fh(a,:,:) = dim_avg_n(track_error_exp(a,:,:,:)-track_error_ctl,1) ; forecast hour only intensity_error_avg_st_fh(a,:,:) = dim_avg_n(intensity_error_exp(a,:,:,:)-intensity_error_ctl,1) ; forecast hour only wind_error_avg_st_fh(a,:,:) = dim_avg_n(wind_error_exp(a,:,:,:)-wind_error_ctl,1) ; forecast hour only, track_error_avg_st_fh(a,:,:) = abs(track_error_exp_st_fh(:,:,a)-track_error_ctl_st_fh) intensity_error_avg_st_fh(a,:,:) = abs(intensity_error_exp_st_fh(:,:,a)-intensity_error_ctl_st_fh) wind_error_avg_st_fh(a,:,:) = abs(wind_error_exp_st_fh(:,:,a)-wind_error_ctl_st_fh) 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 ; check to see if there is missing data, meaning no data for this storm if (.not.all(ismissing(track_error_ctl_st_fh(i,:))) .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,:))).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,:),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,:),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,:),ndtooned(wind_error_exp_st_fh(i,:,:)),0) max_wind = max(wind_temp) min_wind = min(wind_temp) ; determine number of samples do k=0,28 ; loop over forecast hours (output every 6 hours out to 168 hours) fcst_hr_array_ctl2(k) = num(.not.ismissing(track_error_ctl(i,:,k))) end do 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(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,:))) 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(y).le.fcst_hr_array_exp2(a,y)) then number_sample(a,y) = fcst_hr_array_ctl2(y) else number_sample(a,y) = fcst_hr_array_exp2(a,y) 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) = 2.0*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) = 2.0*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) = 2.0*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@gsnScale = True res@trXMinF = 0 res@trXMaxF = x_axis_max 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@tiMainFontHeightF = 0.021 res@xyLineColor = colors_plot(0) plot0 = gsn_csm_xy(wks,x_axis,track_error_ctl_st_fh(i,0:x_axis_max),res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 8.0 plot1 = gsn_csm_xy(wks,x_axis,track_error_exp_st_fh(i,0:x_axis_max,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,track_error_exp_st_fh(i,0:x_axis_max,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,track_error_exp_st_fh(i,0:x_axis_max,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,track_error_exp_st_fh(i,0:x_axis_max,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@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@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 = "# Samples" bres@tiYAxisSide = "Right" bres@tiYAxisAngleF = 270 bres@tmXTOn = False plot5 = gsn_csm_xy(wks,x_axis,fcst_hr_array_ctl2(0:x_axis_max),bres) bres@gsnXYBarChartColors = "firebrick1" 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 = "orange" 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 = "limegreen" 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 = "lightblue" 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 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) end do ; attach significance marker_value = (/3,4,6,9/) do k=0,dimsizes(EXP_name)-1 do w=0,x_axis_max ; loop over forecast hours mres = True 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. track_error_avg_st_fh(k,i,w).gt.interval_track_st_fh(k,i,w)) then ; 95% confidence interval gsn_polymarker(wks,plot0,x_axis(w),5.0,mres) end if end do end do draw(plot0) frame(wks) delete(res) delete(bres) ;;;;;;;;;;;;;;;;;;; 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@gsnScale = True res@trXMinF = 0 res@trXMaxF = x_axis_max 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@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_st_fh(i,0:x_axis_max),res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 8.0 plot1 = gsn_csm_xy(wks,x_axis,intensity_error_exp_st_fh(i,0:x_axis_max,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,intensity_error_exp_st_fh(i,0:x_axis_max,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,intensity_error_exp_st_fh(i,0:x_axis_max,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,intensity_error_exp_st_fh(i,0:x_axis_max,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@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@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 = "# Samples" bres@tiYAxisSide = "Right" bres@tiYAxisAngleF = 270 bres@tmXTOn = False plot5 = gsn_csm_xy(wks,x_axis,fcst_hr_array_ctl2(0:x_axis_max),bres) bres@gsnXYBarChartColors = "firebrick1" 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 = "orange" 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 = "limegreen" 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 = "lightblue" 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 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) end do ; attach significance marker_value = (/3,4,6,9/) do k=0,dimsizes(EXP_name)-1 do w=0,x_axis_max ; loop over forecast hours mres = True 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. intensity_error_avg_st_fh(k,i,w).gt.interval_intensity_st_fh(k,i,w)) then ; 95% confidence interval gsn_polymarker(wks,plot0,x_axis(w),5.0,mres) end if end do end do draw(plot0) frame(wks) delete(res) delete(bres) ; 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@gsnScale = True res@trXMinF = 0 res@trXMaxF = x_axis_max 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 = "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@xyLineThicknessF = 8.0 res@tiMainFontHeightF = 0.021 res@xyLineColor = colors_plot(0) plot0 = gsn_csm_xy(wks,x_axis,wind_error_ctl_st_fh(i,0:x_axis_max),res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 8.0 plot1 = gsn_csm_xy(wks,x_axis,wind_error_exp_st_fh(i,0:x_axis_max,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,wind_error_exp_st_fh(i,0:x_axis_max,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,wind_error_exp_st_fh(i,0:x_axis_max,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,wind_error_exp_st_fh(i,0:x_axis_max,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@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@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 = "# Samples" bres@tiYAxisSide = "Right" bres@tiYAxisAngleF = 270 bres@tmXTOn = False plot5 = gsn_csm_xy(wks,x_axis,fcst_hr_array_ctl2(0:x_axis_max),bres) bres@gsnXYBarChartColors = "firebrick1" 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 = "orange" 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 = "limegreen" 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 = "lightblue" 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 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) end do ; attach significance marker_value = (/3,4,6,9/) do k=0,dimsizes(EXP_name)-1 do w=0,x_axis_max ; loop over forecast hours mres = True 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. wind_error_avg_st_fh(k,i,w).gt.interval_wind_st_fh(k,i,w)) then ; 95% confidence interval gsn_polymarker(wks,plot0,x_axis(w),5.0,mres) end if end do end do draw(plot0) frame(wks) delete(res) delete(bres) 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 end do print("-------------------------------------------------------") print("Total storms analyzed: " + count_storms) print("-------Track Program completed-------------------------") end