load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl" ; This procedure adds metadata (lat, lon, ensemble) to a 3 dimensional variable. undef("createVarMeta") procedure createVarMeta(x:numeric, lat[*]:numeric, lon[*]:numeric, ensemble[*]:numeric) ; add meta data begin x@_FillValue = 1e20 dimx = dimsizes(x) rankx = dimsizes(dimx) if (rankx.eq.3) then x!0 = "lat_0" x!1 = "lon_0" x!2 = "ensemble0" x&lat_0 = lat x&lon_0 = lon x&ensemble0 = ensemble else print("createVarMeta: error: rankx="+rankx) exit end if end ; This function determines how many members from each ensemble system are in each cluster. undef("clusterMembershipCount") function clusterMembershipCount(cluster_membership) begin cluster_mem = where(cluster_membership.ne.1.0,0.0,1.0) cmce_membership = 0 gefs_membership = tointeger(sum(cluster_mem(0:19))) ens_membership = tointeger(sum(cluster_mem(20:69))) return([/cmce_membership,gefs_membership,ens_membership/]) end ; This fuction calculates and plots the cluster forecasts undef("calculate_and_plot_cluster_forecasts") procedure calculate_and_plot_cluster_forecasts(eof_ts, var, var_mean, lat_0, lon_0, ensemble0,var_name,cluster_image_name) begin E0F1_pos_cluster_members = where(eof_ts(0,:).ge.1.0,1.0,eof_ts@_FillValue) EOF1_pos_cluster_members_conform = conform(var, E0F1_pos_cluster_members, (/2/)) EOF1_pos_cluster = where(EOF1_pos_cluster_members_conform .eq. 1.0, var, var@_FillValue) createVarMeta(EOF1_pos_cluster, lat_0, lon_0, ensemble0) EOF1_pos_cluster_mean = dim_avg_Wrap(EOF1_pos_cluster(lat_0|:,lon_0|:,ensemble0|:)) EOF1_pos_cluster_diff = EOF1_pos_cluster_mean - var_mean copy_VarMeta(EOF1_pos_cluster_mean, EOF1_pos_cluster_diff) E0F1_neg_cluster_members = where(eof_ts(0,:).le.-1.0,1.0,eof_ts@_FillValue) EOF1_neg_cluster_members_conform = conform(var, E0F1_neg_cluster_members, (/2/)) EOF1_neg_cluster = where(EOF1_neg_cluster_members_conform .eq. 1.0, var, var@_FillValue) createVarMeta(EOF1_neg_cluster, lat_0, lon_0, ensemble0) EOF1_neg_cluster_mean = dim_avg_Wrap(EOF1_neg_cluster(lat_0|:,lon_0|:,ensemble0|:)) EOF1_neg_cluster_diff = EOF1_neg_cluster_mean - var_mean copy_VarMeta(EOF1_neg_cluster_mean, EOF1_neg_cluster_diff) E0F2_pos_cluster_members = where(eof_ts(1,:).ge.1.0,1.0,eof_ts@_FillValue) EOF2_pos_cluster_members_conform = conform(var, E0F2_pos_cluster_members, (/2/)) EOF2_pos_cluster = where(EOF2_pos_cluster_members_conform .eq. 1.0, var, var@_FillValue) createVarMeta(EOF2_pos_cluster, lat_0, lon_0, ensemble0) EOF2_pos_cluster_mean = dim_avg_Wrap(EOF2_pos_cluster(lat_0|:,lon_0|:,ensemble0|:)) EOF2_pos_cluster_diff = EOF2_pos_cluster_mean - var_mean copy_VarMeta(EOF2_pos_cluster_mean, EOF2_pos_cluster_diff) E0F2_neg_cluster_members = where(eof_ts(1,:).le.-1.0,1.0,eof_ts@_FillValue) EOF2_neg_cluster_members_conform = conform(var, E0F2_neg_cluster_members, (/2/)) EOF2_neg_cluster = where(EOF2_neg_cluster_members_conform .eq. 1.0, var, var@_FillValue) createVarMeta(EOF2_neg_cluster, lat_0, lon_0, ensemble0) EOF2_neg_cluster_mean = dim_avg_Wrap(EOF2_neg_cluster(lat_0|:,lon_0|:,ensemble0|:)) EOF2_neg_cluster_diff = EOF2_neg_cluster_mean - var_mean copy_VarMeta(EOF2_neg_cluster_mean, EOF2_neg_cluster_diff) mean_cluster_members = where((eof_ts(0,:).gt.-1.0) .and. (eof_ts(0,:).lt.1.0) .and. (eof_ts(1,:).gt.-1.0) .and. (eof_ts(1,:).lt.1.0),1.0,eof_ts@_FillValue) mean_cluster_members_conform = conform(var, mean_cluster_members, (/2/)) mean_cluster = where(mean_cluster_members_conform .eq. 1.0, var, var@_FillValue) createVarMeta(mean_cluster, lat_0, lon_0, ensemble0) mean_cluster_mean = dim_avg_Wrap(mean_cluster(lat_0|:,lon_0|:,ensemble0|:)) mean_cluster_diff = mean_cluster_mean - var_mean copy_VarMeta(mean_cluster_mean, mean_cluster_diff) imagedir = "/cf_hmt/hmt/lamberson/day_8_10/fuzzy_clustering/images/cluster_forecast/" cluster_wks_name = imagedir+cluster_image_name cluster_wks = gsn_open_wks("png",cluster_wks_name) gsn_define_colormap(cluster_wks,"amwg_blueyellowred") cluster_diff_res = True cluster_diff_res@gsnDraw = False cluster_diff_res@gsnFrame = False cluster_diff_res@mpProjection = "CylindricalEquidistant" ;cluster_diff_res@mpProjection = "LambertConformal" ;cluster_diff_res@gsnMaskLambertConformal = True cluster_diff_res@mpDataBaseVersion = "MediumRes" cluster_diff_res@mpLimitMode = "LatLon" cluster_diff_res@mpMinLatF = 52. cluster_diff_res@mpMaxLatF = 72. cluster_diff_res@mpMaxLonF = 231. cluster_diff_res@mpMinLonF = 190. cluster_diff_res@gsnAddCyclic = False cluster_diff_res@mpCenterLonF = (cluster_diff_res@mpMaxLonF +cluster_diff_res@mpMinLonF)/2.0 cluster_diff_res@mpLabelsOn = False cluster_diff_res@mpPerimOn = True cluster_diff_res@mpOutlineBoundarySets = "National" cluster_diff_res@mpOutlineSpecifiers = (/"Canada : Provinces","United States : States"/) cluster_diff_res@mpFillOn = False cluster_diff_res@mpGeophysicalLineColor = "saddlebrown" cluster_diff_res@mpGeophysicalLineThicknessF = 3.5 cluster_diff_res@mpNationalLineColor = "saddlebrown" cluster_diff_res@mpNationalLineThicknessF = 3.5 cluster_diff_res@mpUSStateLineColor = "saddlebrown" cluster_diff_res@mpUSStateLineThicknessF = 3.5 cluster_diff_res@gsnMaximize = True cluster_diff_res@gsnPaperOrientation = "portrait" cluster_diff_res@cnLevelSelectionMode = "ExplicitLevels" cluster_diff_res@cnLevels = (/-160.,-140.,-120.,-100.,-80.,-60.,-40.,-20.,20.,40.,60.,80.,100.,120.,140.,160./) cluster_diff_res@cnFillOn = True cluster_diff_res@cnFillColors = (/2,3,4,5,6,7,8,9,0,10,11,12,13,14,15,16,17/) cluster_diff_res@cnLinesOn = False cluster_diff_res@cnInfoLabelOn = False cluster_diff_res@cnLineLabelsOn = False cluster_diff_res@cnInfoLabelOn = False cluster_diff_res@cnLineLabelsOn = False cluster_diff_res@lbLabelBarOn = False cluster_diff_res@gsnLeftString = "" cluster_diff_res@gsnRightString = "" cluster_diff_res@gsnCenterString = "" cluster_diff_res@gsnCenterStringFontHeightF = 0.03 cluster_diff_res@gsnLeftStringFontHeightF = 0.025 cluster_diff_res@gsnRightStringFontHeightF = 0.025 cluster_diff_res@tiMainString = "" ;cluster_diff_res@cnLineDrawOrder = "PostDraw" cluster_mean_res = True cluster_mean_res@gsnDraw = False cluster_mean_res@gsnFrame = False cluster_mean_res@gsnMaximize = True cluster_mean_res@gsnPaperOrientation = "landscape" cluster_mean_res@cnLevelSelectionMode = "ManualLevels" cluster_mean_res@cnMinLevelValF = 4980 cluster_mean_res@cnMaxLevelValF = 5940 cluster_mean_res@cnLevelSpacingF = 60 cluster_mean_res@cnLineColor = "black" cluster_mean_res@cnLineThicknessF = 4 cluster_mean_res@cnLineLabelInterval = 1 cluster_mean_res@cnLineDrawOrder = "PostDraw" cluster_mean_res@cnLineLabelPlacementMode = "Constant" cluster_mean_res@cnLineLabelFont = 22 cluster_mean_res@cnLineLabelFontHeightF = 0.018 cluster_mean_res@gsnLeftString = "" cluster_mean_res@gsnRightString = "" cluster_mean_res@gsnCenterString = "" cluster_mean_res@tiMainString = "" cluster_mean_res@cnInfoLabelOn = False cluster_mean_res@gsnAddCyclic = False cluster_mean_res@gsnContourZeroLineThicknessF = 0.0 cluster_mean_res@gsnContourNegLineDashPattern = 1 if(var_name .eq. "500 hPa Heights") then cluster_diff_res@mpProjection = "CylindricalEquidistant" cluster_diff_res@gsnMaskLambertConformal = False cluster_diff_res@mpMinLatF = 25. cluster_diff_res@mpMaxLatF = 80. cluster_diff_res@mpMaxLonF = 310. cluster_diff_res@mpMinLonF = 190. cluster_mean_res@cnLineLabelInterval = 3 cluster_mean_res@cnLineLabelFontHeightF = 0.012 end if if (var_name .eq. "maximum temperature" .or. var_name .eq. "minimum temperature") then cluster_diff_res@cnLevels = (/-16.,-14.,-12.,-10.,-8.,-6.,-4.,-2.,2.,4.,6.,8.,10.,12.,14.,16./) cluster_mean_res@cnMinLevelValF = -60 cluster_mean_res@cnMaxLevelValF = 110 cluster_mean_res@cnLevelSpacingF = 10 end if if (var_name .eq. "qpf") then cluster_diff_res@cnLevels = (/-2.00,-1.00,-0.75,-0.50,-0.25,-0.10,-0.05,-0.01,0.01,0.05,0.10,0.25,0.50,0.75,1.00,2.00/) cluster_mean_res@cnLevelSelectionMode = "ExplicitLevels" cluster_mean_res@cnLevels = (/0.1,0.25,0.50,0.75,1.0,1.5,2.0,2.5,3.0,4.0/) end if membership = clusterMembershipCount(E0F1_pos_cluster_members) cluster_diff_res@gsnLeftString = "+ EOF1 Cluster" cluster_diff_res@gsnRightString = "C: "+tostring(membership[0])+" G: "+tostring(membership[1])+" E: "+tostring(membership[2]) cluster_diff_res@gsnCenterString = "" EOF1_pos_cluster_diff_plot = gsn_csm_contour_map(cluster_wks,EOF1_pos_cluster_diff,cluster_diff_res) EOF1_pos_cluster_mean_plot = gsn_csm_contour(cluster_wks,EOF1_pos_cluster_mean,cluster_mean_res) overlay(EOF1_pos_cluster_diff_plot,EOF1_pos_cluster_mean_plot) membership = clusterMembershipCount(E0F1_neg_cluster_members) cluster_diff_res@gsnLeftString = "- EOF1 Cluster" cluster_diff_res@gsnRightString = "C: "+tostring(membership[0])+" G: "+tostring(membership[1])+" E: "+tostring(membership[2]) cluster_diff_res@gsnCenterString = "" EOF1_neg_cluster_diff_plot = gsn_csm_contour_map(cluster_wks,EOF1_neg_cluster_diff,cluster_diff_res) EOF1_neg_cluster_mean_plot = gsn_csm_contour(cluster_wks,EOF1_neg_cluster_mean,cluster_mean_res) overlay(EOF1_neg_cluster_diff_plot,EOF1_neg_cluster_mean_plot) membership = clusterMembershipCount(mean_cluster_members) cluster_diff_res@gsnLeftString = "Mean Cluster" cluster_diff_res@gsnRightString = "C: "+tostring(membership[0])+" G: "+tostring(membership[1])+" E: "+tostring(membership[2]) cluster_diff_res@gsnCenterString = "" MEAN_cluster_diff_plot = gsn_csm_contour_map(cluster_wks,mean_cluster_diff,cluster_diff_res) MEAN_cluster_mean_plot = gsn_csm_contour(cluster_wks,mean_cluster_mean,cluster_mean_res) overlay(MEAN_cluster_diff_plot,MEAN_cluster_mean_plot) membership = clusterMembershipCount(E0F2_pos_cluster_members) cluster_diff_res@gsnLeftString = "+ EOF2 Cluster" cluster_diff_res@gsnRightString = "C: "+tostring(membership[0])+" G: "+tostring(membership[1])+" E: "+tostring(membership[2]) cluster_diff_res@gsnCenterString = "" EOF2_pos_cluster_diff_plot = gsn_csm_contour_map(cluster_wks,EOF2_pos_cluster_diff,cluster_diff_res) EOF2_pos_cluster_mean_plot = gsn_csm_contour(cluster_wks,EOF2_pos_cluster_mean,cluster_mean_res) overlay(EOF2_pos_cluster_diff_plot,EOF2_pos_cluster_mean_plot) membership = clusterMembershipCount(E0F2_neg_cluster_members) cluster_diff_res@gsnLeftString = "- EOF2 Cluster" cluster_diff_res@gsnRightString = "C: "+tostring(membership[0])+" G: "+tostring(membership[1])+" E: "+tostring(membership[2]) cluster_diff_res@gsnCenterString = "" EOF2_neg_cluster_diff_plot = gsn_csm_contour_map(cluster_wks,EOF2_neg_cluster_diff,cluster_diff_res) EOF2_neg_cluster_mean_plot = gsn_csm_contour(cluster_wks,EOF2_neg_cluster_mean,cluster_mean_res) overlay(EOF2_neg_cluster_diff_plot,EOF2_neg_cluster_mean_plot) cluster_hgt_plots = new(5,graphic) cluster_hgt_plots(0) = EOF1_pos_cluster_diff_plot cluster_hgt_plots(1) = EOF1_neg_cluster_diff_plot cluster_hgt_plots(2) = MEAN_cluster_diff_plot cluster_hgt_plots(3) = EOF2_pos_cluster_diff_plot cluster_hgt_plots(4) = EOF2_neg_cluster_diff_plot panel_res = True panel_res@gsnMaximize = True ;panel_res@gsnFrame = False panel_res@gsnPanelRowSpec = True panel_res@gsnPanelLabelBar = True panel_res@lbLabelFontHeightF = 0.011 panel_res@lbLabelStride = 1 panel_res@gsnPaperOrientation = "landscape" panel_res@pmLabelBarHeightF = 0.08 panel_res@pmLabelBarWidthF = 0.70 gsn_panel(cluster_wks, cluster_hgt_plots,(/2,1,2/),panel_res) delete(cluster_wks) delete(cluster_hgt_plots) end ; This procedure calculates and plots the cluster forecasts of maximum temperature undef("max_temp_cluster_forecasts") procedure max_temp_cluster_forecasts(eof_ts) begin model_dir = "/cf_hmt/hmt/lamberson/day_8_10/fuzzy_clustering/model_data/cluster/" ;cmce_files = systemfunc("ls "+model_dir+"cmc_ge*.grb") ;cmce_f = addfiles(cmce_files, "r") ;ListSetType(cmce_f, "join") ;cmce_var_3d = cmce_f[:]->TMAX_P11_L103_GLL0_max6h(:,::-1,:) ;cmce_var_4d = reshape(cmce_var_3d,(/20,12,361,720/)) gefs_files = systemfunc("ls "+model_dir+"ge*.grb") gefs_f = addfiles(gefs_files, "r") ListSetType(gefs_f, "join") gefs_var_3d = gefs_f[:]->TMAX_P11_L103_GLL0_max6h gefs_var_4d = reshape(gefs_var_3d,(/20,12,361,720/)) ens_files = systemfunc("ls "+model_dir+"E1E*.grb") ens_f = addfiles(ens_files, "r") ListSetType(ens_f, "join") ens_var_4d_ugh = ens_f[1:12]->MX2T6_GDS0_SFC(:,1:50,:,:) ens_var_4d = ens_var_4d_ugh ens_var_4d(:,:,:,0:359) = ens_var_4d_ugh(:,:,:,360:719) ens_var_4d(:,:,:,360:719) = ens_var_4d_ugh(:,:,:,0:359) lat_0 = gefs_f[0]->lat_0 lon_0 = gefs_f[0]->lon_0 ensemble0 = ispan(0,69,1) j = 0 do i = 0,2 ;cmce_max_temps = dim_max_n(cmce_var_4d(:,j:j+3,:,:),(/1/)) gefs_max_temps = dim_max_n(gefs_var_4d(:,j:j+3,:,:),(/1/)) ens_max_temps = dim_max_n(ens_var_4d(j:j+3,:,:,:),(/0/)) max_temp_combined_1 = array_append_record(gefs_max_temps,ens_max_temps,0) printVarSummary(max_temp_combined_1) max_temp = max_temp_combined_1(dim1|:,dim2|:,dim0|:) max_temp = 1.8* (max_temp - 273.15) + 32 createVarMeta(max_temp, lat_0, lon_0, ensemble0) max_temp_mean = dim_avg_Wrap(max_temp(lat_0|:,lon_0|:,ensemble0|:)) if(i .eq. 0) then cluster_image_name = "TMAX/day_8_cluster_maximum_temperatures_ak.png" else if(i .eq. 1) then cluster_image_name = "TMAX/day_9_cluster_maximum_temperatures_ak.png" else cluster_image_name = "TMAX/day_10_cluster_maximum_temperatures_ak.png" end if end if calculate_and_plot_cluster_forecasts(eof_ts,max_temp,max_temp_mean,lat_0,lon_0,ensemble0,"maximum temperature",cluster_image_name) j = j+4 delete(max_temp) end do end undef("min_temp_cluster_forecasts") procedure min_temp_cluster_forecasts(eof_ts) begin model_dir = "/cf_hmt/hmt/lamberson/day_8_10/fuzzy_clustering/model_data/cluster/" ;cmce_files = systemfunc("ls "+model_dir+"cmc_ge*.grb") ;cmce_f = addfiles(cmce_files, "r") ;ListSetType(cmce_f, "join") ;cmce_var_3d = cmce_f[:]->TMIN_P11_L103_GLL0_min6h(:,::-1,:) ;cmce_var_4d = reshape(cmce_var_3d,(/20,12,361,720/)) gefs_files = systemfunc("ls "+model_dir+"ge*.grb") gefs_f = addfiles(gefs_files, "r") ListSetType(gefs_f, "join") gefs_var_3d = gefs_f[:]->TMIN_P11_L103_GLL0_min6h gefs_var_4d = reshape(gefs_var_3d,(/20,12,361,720/)) ens_files = systemfunc("ls "+model_dir+"E1E*.grb") ens_f = addfiles(ens_files, "r") ListSetType(ens_f, "join") ens_var_4d_ugh = ens_f[1:12]->MN2T6_GDS0_SFC(:,1:50,:,:) ens_var_4d = ens_var_4d_ugh ens_var_4d(:,:,:,0:359) = ens_var_4d_ugh(:,:,:,360:719) ens_var_4d(:,:,:,360:719) = ens_var_4d_ugh(:,:,:,0:359) lat_0 = gefs_f[0]->lat_0 lon_0 = gefs_f[0]->lon_0 ensemble0 = ispan(0,69,1) j = 0 do i = 0,2 ;cmce_min_temps = dim_min_n(cmce_var_4d(:,j:j+3,:,:),(/1/)) gefs_min_temps = dim_min_n(gefs_var_4d(:,j:j+3,:,:),(/1/)) ens_min_temps = dim_min_n(ens_var_4d(j:j+3,:,:,:),(/0/)) min_temp_combined_1 = array_append_record(gefs_min_temps,ens_min_temps,0) ;min_temp_combined_2 = array_append_record(min_temp_combined_1,ens_min_temps,0) min_temp = min_temp_combined_1(dim1|:,dim2|:,dim0|:) min_temp = 1.8* (min_temp - 273.15) + 32 createVarMeta(min_temp, lat_0, lon_0, ensemble0) min_temp_mean = dim_avg_Wrap(min_temp(lat_0|:,lon_0|:,ensemble0|:)) if(i .eq. 0) then cluster_image_name = "TMIN/day_8_cluster_minimum_temperatures_ak.png" else if(i .eq. 1) then cluster_image_name = "TMIN/day_9_cluster_minimum_temperatures_ak.png" else cluster_image_name = "TMIN/day_10_cluster_minimum_temperatures_ak.png" end if end if calculate_and_plot_cluster_forecasts(eof_ts,min_temp,min_temp_mean,lat_0,lon_0,ensemble0,"minimum temperature",cluster_image_name) j = j+4 delete(min_temp) end do end undef("qpf_cluster_forecasts") procedure qpf_cluster_forecasts(eof_ts) begin model_dir = "/cf_hmt/hmt/lamberson/day_8_10/fuzzy_clustering/model_data/cluster/" ;cmce_files = systemfunc("ls "+model_dir+"cmc_ge*.grb") ;cmce_f = addfiles(cmce_files, "r") ;ListSetType(cmce_f, "join") ;cmce_var_3d = cmce_f[:]->APCP_P11_L1_GLL0_acc6h(:,::-1,:) ;cmce_var_4d = reshape(cmce_var_3d,(/20,12,361,720/)) gefs_files = systemfunc("ls "+model_dir+"ge*.grb") gefs_f = addfiles(gefs_files, "r") ListSetType(gefs_f, "join") gefs_var_3d = gefs_f[:]->APCP_P11_L1_GLL0_acc6h gefs_var_4d = reshape(gefs_var_3d,(/20,12,361,720/)) ens_files = systemfunc("ls "+model_dir+"E1E*.grb") ens_f = addfiles(ens_files, "r") ListSetType(ens_f, "join") ens_var_4d_ugh = ens_f[:]->TP_GDS0_SFC(:,1:50,:,:) ens_var_4d = ens_var_4d_ugh ens_var_4d(:,:,:,0:359) = ens_var_4d_ugh(:,:,:,360:719) ens_var_4d(:,:,:,360:719) = ens_var_4d_ugh(:,:,:,0:359) lat_0 = gefs_f[0]->lat_0 lon_0 = gefs_f[0]->lon_0 ensemble0 = ispan(0,69,1) ;cmce_qpf = dim_sum_n(cmce_var_4d,(/1/)) gefs_qpf = dim_sum_n(gefs_var_4d,(/1/)) ens_qpf = ens_var_4d(12,:,:,:) - ens_var_4d(0,:,:,:) ;cmce_qpf = cmce_qpf*(1.0/25.4) gefs_qpf = gefs_qpf*(1.0/25.4) ens_qpf = ens_qpf*(1.0/0.0254) qpf_combined_1 = array_append_record(gefs_qpf,ens_qpf,0) ;qpf_combined_2 = array_append_record(qpf_combined_1,ens_qpf,0) qpf = qpf_combined_1(dim1|:,dim2|:,dim0|:) createVarMeta(qpf, lat_0, lon_0, ensemble0) qpf_mean = dim_avg_Wrap(qpf(lat_0|:,lon_0|:,ensemble0|:)) cluster_image_name = "QPF/72_hour_cluster_qpf_ak.png" calculate_and_plot_cluster_forecasts(eof_ts,qpf,qpf_mean,lat_0,lon_0,ensemble0,"qpf",cluster_image_name) delete(qpf) end ; This fuction calculates and plots the evolution of each cluster from forecast hour 0 undef("calculate_and_plot_cluster_evolution") procedure calculate_and_plot_cluster_evolution(eof_ts) begin model_dir = "/cf_hmt/hmt/lamberson/day_8_10/fuzzy_clustering/model_data/cluster_evolution/" cmce_files_evolution = systemfunc("ls "+model_dir+"cmc_ge*.grb") print(cmce_files_evolution) cmce_f_evolution = addfiles(cmce_files_evolution, "r") ListSetType(cmce_f_evolution, "join") cmce_var_3d_evolution = cmce_f_evolution[:]->HGT_P1_L100_GLL0(:,5,::-1,:) cmce_var_4d_evolution = reshape(cmce_var_3d_evolution,(/20,16,361,720/)) cmce_var_file_evolution_4d = reshape(cmce_files_evolution,(/20,16/)) print(cmce_var_file_evolution_4d) printVarSummary(cmce_var_4d_evolution) gefs_files_evolution = systemfunc("ls "+model_dir+"ge*.grb") gefs_f_evolution = addfiles(gefs_files_evolution, "r") ListSetType(gefs_f_evolution, "join") gefs_var_3d_evolution = gefs_f_evolution[:]->HGT_P1_L100_GLL0(:,5,:,:) gefs_var_4d_evolution = reshape(gefs_var_3d_evolution,(/20,16,361,720/)) ens_files_evolution = systemfunc("ls "+model_dir+"E1E*.grb") ens_f_evolution = addfiles(ens_files_evolution, "r") ListSetType(ens_f_evolution, "join") ens_var_4d_evolution_temp = ens_f_evolution[:]->GH_GDS0_ISBL(:,1:50,2,:,:) ens_var_4d_evolution = ens_var_4d_evolution_temp(ensemble0|:,ncl_join|:,g0_lat_2|:,g0_lon_3|:) printVarSummary(ens_var_4d_evolution) lat_0 = gefs_f_evolution[0]->lat_0 lon_0 = gefs_f_evolution[0]->lon_0 ensemble0 = ispan(0,89,1) var_evolution_combined_1 = array_append_record(cmce_var_4d_evolution,gefs_var_4d_evolution,0) var_evolution_combined_2 = array_append_record(var_evolution_combined_1,ens_var_4d_evolution,0) printVarSummary(var_evolution_combined_2) do i = 0,15 var_evolution_3d = var_evolution_combined_2(:,i,:,:) var_evolution = var_evolution_3d(dim2|:,dim3|:,dim0|:) printVarSummary(var_evolution) var_evolution_mean = dim_avg_n_Wrap(var_evolution,(/2/)) E0F1_pos_cluster_members = where(eof_ts(0,:).ge.1.0,1.0,eof_ts@_FillValue) EOF1_pos_cluster_members_conform = conform(var_evolution, E0F1_pos_cluster_members, (/2/)) EOF1_pos_cluster = where(EOF1_pos_cluster_members_conform .eq. 1.0, var_evolution, var_evolution@_FillValue) createVarMeta(EOF1_pos_cluster, lat_0, lon_0, ensemble0) EOF1_pos_cluster_mean = dim_avg_Wrap(EOF1_pos_cluster(lat_0|:,lon_0|:,ensemble0|:)) EOF1_pos_cluster_diff = EOF1_pos_cluster_mean - var_evolution_mean copy_VarMeta(EOF1_pos_cluster_mean, EOF1_pos_cluster_diff) E0F1_neg_cluster_members = where(eof_ts(0,:).le.-1.0,1.0,eof_ts@_FillValue) EOF1_neg_cluster_members_conform = conform(var_evolution, E0F1_neg_cluster_members, (/2/)) EOF1_neg_cluster = where(EOF1_neg_cluster_members_conform .eq. 1.0, var_evolution, var_evolution@_FillValue) createVarMeta(EOF1_neg_cluster, lat_0, lon_0, ensemble0) EOF1_neg_cluster_mean = dim_avg_Wrap(EOF1_neg_cluster(lat_0|:,lon_0|:,ensemble0|:)) EOF1_neg_cluster_diff = EOF1_neg_cluster_mean - var_evolution_mean copy_VarMeta(EOF1_neg_cluster_mean, EOF1_neg_cluster_diff) E0F2_pos_cluster_members = where(eof_ts(1,:).ge.1.0,1.0,eof_ts@_FillValue) EOF2_pos_cluster_members_conform = conform(var_evolution, E0F2_pos_cluster_members, (/2/)) EOF2_pos_cluster = where(EOF2_pos_cluster_members_conform .eq. 1.0, var_evolution, var_evolution@_FillValue) createVarMeta(EOF2_pos_cluster, lat_0, lon_0, ensemble0) EOF2_pos_cluster_mean = dim_avg_Wrap(EOF2_pos_cluster(lat_0|:,lon_0|:,ensemble0|:)) EOF2_pos_cluster_diff = EOF2_pos_cluster_mean - var_evolution_mean copy_VarMeta(EOF2_pos_cluster_mean, EOF2_pos_cluster_diff) E0F2_neg_cluster_members = where(eof_ts(1,:).le.-1.0,1.0,eof_ts@_FillValue) EOF2_neg_cluster_members_conform = conform(var_evolution, E0F2_neg_cluster_members, (/2/)) EOF2_neg_cluster = where(EOF2_neg_cluster_members_conform .eq. 1.0, var_evolution, var_evolution@_FillValue) createVarMeta(EOF2_neg_cluster, lat_0, lon_0, ensemble0) EOF2_neg_cluster_mean = dim_avg_Wrap(EOF2_neg_cluster(lat_0|:,lon_0|:,ensemble0|:)) EOF2_neg_cluster_diff = EOF2_neg_cluster_mean - var_evolution_mean copy_VarMeta(EOF2_neg_cluster_mean, EOF2_neg_cluster_diff) mean_cluster_members = where((eof_ts(0,:).gt.-1.0) .and. (eof_ts(0,:).lt.1.0) .and. (eof_ts(1,:).gt.-1.0) .and. (eof_ts(1,:).lt.1.0),1.0,eof_ts@_FillValue) mean_cluster_members_conform = conform(var_evolution, mean_cluster_members, (/2/)) mean_cluster = where(mean_cluster_members_conform .eq. 1.0, var_evolution, var_evolution@_FillValue) createVarMeta(mean_cluster, lat_0, lon_0, ensemble0) mean_cluster_mean = dim_avg_Wrap(mean_cluster(lat_0|:,lon_0|:,ensemble0|:)) mean_cluster_diff = mean_cluster_mean - var_evolution_mean copy_VarMeta(mean_cluster_mean, mean_cluster_diff) cluster_diff_res = True cluster_diff_res@gsnDraw = False cluster_diff_res@gsnFrame = False cluster_diff_res@mpProjection = "CylindricalEquidistant" cluster_diff_res@mpDataBaseVersion = "MediumRes" cluster_diff_res@mpLimitMode = "LatLon" cluster_diff_res@mpMinLatF = 10. cluster_diff_res@mpMaxLatF = 90. cluster_diff_res@mpMaxLonF = 330. cluster_diff_res@mpMinLonF = 100. cluster_diff_res@gsnAddCyclic = False cluster_diff_res@mpCenterLonF = (cluster_diff_res@mpMaxLonF +cluster_diff_res@mpMinLonF)/2.0 cluster_diff_res@mpLabelsOn = False cluster_diff_res@mpPerimOn = True cluster_diff_res@mpOutlineBoundarySets = "National" cluster_diff_res@mpOutlineSpecifiers = (/"Canada : Provinces","United States : States"/) cluster_diff_res@mpFillOn = False cluster_diff_res@mpGeophysicalLineColor = "saddlebrown" cluster_diff_res@mpGeophysicalLineThicknessF = 3.5 cluster_diff_res@mpNationalLineColor = "saddlebrown" cluster_diff_res@mpNationalLineThicknessF = 3.5 cluster_diff_res@mpUSStateLineColor = "saddlebrown" cluster_diff_res@mpUSStateLineThicknessF = 3.5 cluster_diff_res@gsnMaximize = True cluster_diff_res@gsnPaperOrientation = "portrait" cluster_diff_res@cnLevelSelectionMode = "ExplicitLevels" cluster_diff_res@cnLevels = (/-160.,-140.,-120.,-100.,-80.,-60.,-40.,-20.,20.,40.,60.,80.,100.,120.,140.,160./) cluster_diff_res@cnFillOn = True cluster_diff_res@cnFillColors = (/2,3,4,5,6,7,8,9,0,10,11,12,13,14,15,16,17/) cluster_diff_res@cnLinesOn = False cluster_diff_res@cnInfoLabelOn = False cluster_diff_res@cnLineLabelsOn = False cluster_diff_res@cnInfoLabelOn = False cluster_diff_res@cnLineLabelsOn = False cluster_diff_res@lbLabelBarOn = True cluster_diff_res@gsnLeftString = "" cluster_diff_res@gsnRightString = "" cluster_diff_res@gsnCenterString = "+ EOF1 Cluster" cluster_diff_res@gsnCenterStringFontHeightF = 0.03 cluster_diff_res@gsnLeftStringFontHeightF = 0.01 cluster_diff_res@gsnRightStringFontHeightF = 0.01 cluster_diff_res@tiMainString = "" cluster_mean_res = True cluster_mean_res@gsnDraw = False cluster_mean_res@gsnFrame = False cluster_mean_res@gsnMaximize = True cluster_mean_res@gsnPaperOrientation = "landscape" cluster_mean_res@cnLevelSelectionMode = "ManualLevels" cluster_mean_res@cnMinLevelValF = 4980 cluster_mean_res@cnMaxLevelValF = 5940 cluster_mean_res@cnLevelSpacingF = 60 cluster_mean_res@cnLineColor = "black" cluster_mean_res@cnLineThicknessF = 2.2 cluster_mean_res@cnLineLabelInterval = 2 cluster_mean_res@cnLineDrawOrder = "PostDraw" cluster_mean_res@cnLineLabelPlacementMode = "Constant" cluster_mean_res@cnLineLabelFont = 22 cluster_mean_res@cnLineLabelFontHeightF = 0.008 cluster_mean_res@gsnLeftString = "" cluster_mean_res@gsnRightString = "" cluster_mean_res@gsnCenterString = "" cluster_mean_res@tiMainString = "" cluster_mean_res@cnInfoLabelOn = False cluster_mean_res@gsnAddCyclic = False cluster_mean_res@gsnContourZeroLineThicknessF = 2.2 cluster_mean_res@gsnContourNegLineDashPattern = 1 imagedir = "/cf_hmt/hmt/lamberson/day_8_10/fuzzy_clustering/images/cluster_evolution/" cluster_image_name = "pos_EOF1/500_HGT_"+tostring(i) cluster_wks_name = imagedir+cluster_image_name cluster_wks = gsn_open_wks("png",cluster_wks_name) gsn_define_colormap(cluster_wks,"amwg_blueyellowred") EOF1_pos_cluster_diff_plot = gsn_csm_contour_map(cluster_wks,EOF1_pos_cluster_diff,cluster_diff_res) EOF1_pos_cluster_mean_plot = gsn_csm_contour(cluster_wks,EOF1_pos_cluster_mean,cluster_mean_res) overlay(EOF1_pos_cluster_diff_plot,EOF1_pos_cluster_mean_plot) draw(EOF1_pos_cluster_diff_plot) frame(cluster_wks) delete(cluster_wks) cluster_image_name = "neg_EOF1/500_HGT_"+tostring(i) cluster_wks_name = imagedir+cluster_image_name cluster_wks = gsn_open_wks("png",cluster_wks_name) gsn_define_colormap(cluster_wks,"amwg_blueyellowred") cluster_diff_res@gsnCenterString = "- EOF1 Cluster" EOF1_neg_cluster_diff_plot = gsn_csm_contour_map(cluster_wks,EOF1_neg_cluster_diff,cluster_diff_res) EOF1_neg_cluster_mean_plot = gsn_csm_contour(cluster_wks,EOF1_neg_cluster_mean,cluster_mean_res) overlay(EOF1_neg_cluster_diff_plot,EOF1_neg_cluster_mean_plot) draw(EOF1_neg_cluster_diff_plot) frame(cluster_wks) delete(cluster_wks) cluster_image_name = "MEAN/500_HGT_"+tostring(i) cluster_wks_name = imagedir+cluster_image_name cluster_wks = gsn_open_wks("png",cluster_wks_name) gsn_define_colormap(cluster_wks,"amwg_blueyellowred") cluster_diff_res@gsnCenterString = "Mean Cluster" MEAN_cluster_diff_plot = gsn_csm_contour_map(cluster_wks,mean_cluster_diff,cluster_diff_res) MEAN_cluster_mean_plot = gsn_csm_contour(cluster_wks,mean_cluster_mean,cluster_mean_res) overlay(MEAN_cluster_diff_plot,MEAN_cluster_mean_plot) draw(MEAN_cluster_diff_plot) frame(cluster_wks) delete(cluster_wks) cluster_image_name = "pos_EOF2/500_HGT_"+tostring(i) cluster_wks_name = imagedir+cluster_image_name cluster_wks = gsn_open_wks("png",cluster_wks_name) gsn_define_colormap(cluster_wks,"amwg_blueyellowred") cluster_diff_res@gsnCenterString = "+ EOF2 Cluster" EOF2_pos_cluster_diff_plot = gsn_csm_contour_map(cluster_wks,EOF2_pos_cluster_diff,cluster_diff_res) EOF2_pos_cluster_mean_plot = gsn_csm_contour(cluster_wks,EOF2_pos_cluster_mean,cluster_mean_res) overlay(EOF2_pos_cluster_diff_plot,EOF2_pos_cluster_mean_plot) draw(EOF2_pos_cluster_diff_plot) frame(cluster_wks) delete(cluster_wks) cluster_image_name = "neg_EOF2/500_HGT_"+tostring(i) cluster_wks_name = imagedir+cluster_image_name cluster_wks = gsn_open_wks("png",cluster_wks_name) gsn_define_colormap(cluster_wks,"amwg_blueyellowred") cluster_diff_res@gsnCenterString = "- EOF2 Cluster" EOF2_neg_cluster_diff_plot = gsn_csm_contour_map(cluster_wks,EOF2_neg_cluster_diff,cluster_diff_res) EOF2_neg_cluster_mean_plot = gsn_csm_contour(cluster_wks,EOF2_neg_cluster_mean,cluster_mean_res) overlay(EOF2_neg_cluster_diff_plot,EOF2_neg_cluster_mean_plot) draw(EOF2_neg_cluster_diff_plot) frame(cluster_wks) delete(cluster_wks) ; delete(cluster_wks) ; delete(cluster_hgt_plots) end do end begin imagedir = "/cf_hmt/hmt/lamberson/day_8_10/fuzzy_clustering/images/" model_dir = "/cf_hmt/hmt/lamberson/day_8_10/fuzzy_clustering/model_data/cluster/" ;cmce_files = systemfunc("ls "+model_dir+"cmc_ge*.grb") ;cmce_f = addfiles(cmce_files, "r") ;ListSetType(cmce_f, "join") ;cmce_var_3d = cmce_f[:]->HGT_P1_L100_GLL0(:,6,::-1,:) ;cmce_var_4d = reshape(cmce_var_3d,(/20,12,361,720/)) ;cmce_var = dim_avg_n_Wrap(cmce_var_4d,(/1/)) gefs_files = systemfunc("ls "+model_dir+"ge*.grb") gefs_f = addfiles(gefs_files, "r") ListSetType(gefs_f, "join") gefs_var_3d = gefs_f[:]->HGT_P1_L100_GLL0(:,6,:,:) gefs_var_4d = reshape(gefs_var_3d,(/20,12,361,720/)) gefs_var = dim_avg_n_Wrap(gefs_var_4d,(/1/)) ens_files = systemfunc("ls "+model_dir+"E1E*.grb") ens_f = addfiles(ens_files, "r") ListSetType(ens_f, "join") ens_var_4d = ens_f[1:12]->GH_GDS0_ISBL(:,1:50,5,:,:) ens_var_ugh = dim_avg_n_Wrap(ens_var_4d,(/0/)) ens_var = ens_var_ugh ens_var(:,:,0:359) = ens_var_ugh(:,:,360:719) ens_var(:,:,360:719) = ens_var_ugh(:,:,0:359) dummy_combined_var_1 = array_append_record(gefs_var,ens_var,0) printVarSummary(dummy_combined_var_1) ;dummy_combined_var_2 = array_append_record(dummy_combined_var_1,ens_var,0) var = dummy_combined_var_1(g0_lat_1|:,g0_lon_2|:,ensemble0|:) printVarSummary(var) lat_0 = gefs_f[0]->lat_0 lon_0 = gefs_f[0]->lon_0 ensemble0 = ispan(0,69,1) createVarMeta(var, lat_0, lon_0, ensemble0) printVarSummary(var) var_mean = dim_avg_Wrap(var(lat_0|:,lon_0|:,ensemble0|:)) var_mean_conform = conform(var, var_mean,(/0,1/)) printVarSummary(var) ;var = var-273.15 var_normalized = (var-var_mean_conform) wgt = sqrt(cos(lat_0*0.01745329)) var_normalized = var_normalized*conform(var_normalized, wgt, 0) createVarMeta(var_normalized, lat_0, lon_0, ensemble0) ;max_lat = 80.0 ;min_lat = 25.0 ;max_lon = 310.0 ;min_lon = 190.0 max_lat = 80.0 min_lat = 20.0 max_lon = 330.0 min_lon = 180.0 eofmax = 2 eof = eofunc_Wrap(var_normalized({min_lat:max_lat},{min_lon:max_lon},:),eofmax,False) eof_ts = eofunc_ts_Wrap(var_normalized({min_lat:max_lat},{min_lon:max_lon},:),eof,False) do e = 0,eofmax-1 eof(e,:,:) = eof(e,:,:)*sqrt(eof@eval(e)) ; units same as data eof_ts(e,:) = eof_ts(e,:)/sqrt(eof@eval(e)) end do print(eof_ts) wks_name=imagedir+"EOFs.png" wks = gsn_open_wks ("png", wks_name ) gsn_define_colormap(wks, "amwg_blueyellowred") res = True res@gsnFrame = False res@gsnDraw = False res@gsnMaximize = True res@gsnLeftString = "" res@gsnRightString = "" res@gsnCenterString = "" eof_res = res eof_res@mpProjection = "CylindricalEquidistant" eof_res@mpDataBaseVersion = "MediumRes" eof_res@mpLimitMode = "LatLon" eof_res@mpMinLatF = min_lat eof_res@mpMaxLatF = max_lat eof_res@mpMinLonF = min_lon eof_res@mpMaxLonF = max_lon eof_res@gsnAddCyclic = False eof_res@mpCenterLonF = (eof_res@mpMinLonF+eof_res@mpMaxLonF)/2.0 eof_res@mpOutlineBoundarySets = "National" eof_res@mpOutlineSpecifiers = (/"Canada: Provinces","United States : States"/) eof_res@mpFillOn = False eof_res@mpGeophysicalLineColor = "black" eof_res@mpGeophysicalLineThicknessF = 3.0 eof_res@mpNationalLineColor = "black" eof_res@mpNationalLineThicknessF = 3.0 eof_res@mpUSStateLineColor = "black" eof_res@mpUSStateLineThicknessF = 3.0 eof_res@cnLevelSelectionMode = "ExplicitLevels" ;eof_res@cnLevels = (/-0.8,-0.6,-0.4,-0.2,0.2,0.4,0.6,0.8/) eof_res@cnLevels = (/-160,-140,-120,-100,-80,-60,-40,-20,20,40,60,80,100,120,140,160/) ;eof_res@cnLevels = (/-8,-7,-6,-5,-4,-3,-2,-1,1,2,3,4,5,6,7,8/) eof_res@cnFillColors = (/2,3,4,5,6,7,8,9,0,10,11,12,13,14,15,16,17/) eof_res@cnFillOn = True eof_res@cnLinesOn = False eof_res@cnInfoLabelOn = False eof_res@cnLineLabelsOn = False eof_res@lbLabelBarOn = False eof_res@tiMainString = "" eof_res@gsnStringFont = 22 eof_res@gsnLeftString = "EOF1 of Day 8-10 Mean 500 hPa Height Field" eof_res@gsnRightString = "Percent Variance = "+sprintf("%5.1f", eof@pcvar(0)) +"%" eof_res@lbLabelStride = 1 mean_res = res mean_res@cnFillOn = False ;mean_res@cnLevelSelectionMode = "AutomaticLevels" mean_res@cnLevelSelectionMode = "ManualLevels" mean_res@cnMinLevelValF = 4980 mean_res@cnMaxLevelValF = 5940 mean_res@cnLevelSpacingF = 60 ;mean_res@cnMinLevelValF = -40 ;mean_res@cnMaxLevelValF = 40 ;mean_res@cnLevelSpacingF = 2 mean_res@cnLineColor = "grey" mean_res@cnLineThicknessF = 5 mean_res@cnLineLabelInterval = 2 mean_res@cnLineLabelPlacementMode = "Constant" mean_res@cnLineLabelFont = 22 mean_res@cnLineLabelFontHeightF = 0.010 mean_res@cnInfoLabelOn = False EOF1_plot = gsn_csm_contour_map(wks,eof(0,:,:),eof_res) EOF1_plot_ensemble_mean = gsn_csm_contour(wks,var_mean,mean_res) overlay(EOF1_plot,EOF1_plot_ensemble_mean) eof_res@gsnLeftString = "EOF2 of Day 8-10 Mean 500 hPa Height Field" eof_res@gsnRightString = "Percent Variance = "+sprintf("%5.1f", eof@pcvar(1)) +"%" EOF2_plot = gsn_csm_contour_map(wks,eof(1,:,:),eof_res) EOF2_plot_ensemble_mean = gsn_csm_contour(wks,var_mean,mean_res) overlay(EOF2_plot,EOF2_plot_ensemble_mean) plots = new(2,graphic) plots(0) = EOF1_plot plots(1) = EOF2_plot panel_res = True panel_res@gsnMaximize = True ;panel_res@gsnFrame = False panel_res@gsnPanelRowSpec = True panel_res@gsnPanelLabelBar = True panel_res@lbLabelFontHeightF = 0.012 panel_res@lbLabelStride = 1 panel_res@gsnPaperOrientation = "landscape" panel_res@pmLabelBarHeightF = 0.08 panel_res@pmLabelBarWidthF = 0.70 gsn_panel(wks, plots,(/1,1/),panel_res) wks_name_2=imagedir+"cluster_phase_space.png" wks_2 = gsn_open_wks ("png",wks_name_2) phase_res = res phase_res@gsnFrame = False phase_res@gsnDraw = False phase_res@gsnMaximize = True phase_res@tiMainString = "Scatter Plot" phase_res@xyMarkLineMode = "Markers" phase_res@xyMarkers = 5 phase_res@xyMarkerColor = "green" phase_res@xyMarkerSizeF = 0.012 phase_res@xyMarkerThicknessF = 6.0 phase_res@gsnYRefLine = 0.0 phase_res@gsnXRefLine = 0.0 phase_res@trYMinF = -3.5 phase_res@trYMaxF = 3.5 phase_res@trXMinF = -3.5 phase_res@trXMaxF = 3.5 phase_res@tiXAxisString = "PC 1" phase_res@tiYAxisString = "PC 2" phase_res@tmXMajorGrid = True phase_res@tmYMajorGrid = True phase_res@tmXMajorGridLineDashPattern = 2 phase_res@tmYMajorGridLineDashPattern = 2 phase_res@tmXBMode = "Manual" phase_res@tmXBTickStartF = -3.5 phase_res@tmXBTickEndF = 3.5 phase_res@tmXBTickSpacingF = 0.5 phase_res@tmXBLabelFontHeightF = 0.012 phase_res@tmYLMode = "Manual" phase_res@tmYLTickStartF = -3.5 phase_res@tmYLTickEndF = 3.5 phase_res@tmYLTickSpacingF = 0.5 phase_res@tmYLLabelFontHeightF = 0.012 phase_res@tmXTMinorOn = False phase_res@tmXBMinorOn = False phase_res@tmYLMinorOn = False phase_res@tmYRMinorOn = False res_pyline = True res_pyline@gsLineColor = "black" res_pyline@gsLineThicknessF = 10.0 labels1 = (/"GEFS","ECMWF"/) colors1 = (/"green","blue"/) markers1 = (/5, 5/) xpos1 = (/0.50, 0.65/) xpos2 = (/0.53, 0.68/) mkres = True ; Marker resources txres = True ; Text resources txres@txFontHeightF = 0.015 txres@txFont = 22 txres@txJust = "CenterLeft" do i = 0,1 mkres@gsMarkerThicknessF = 6.0 mkres@gsMarkerSizeF = 0.015 mkres@gsMarkerIndex = markers1(i) mkres@gsMarkerColor = colors1(i) gsn_polymarker_ndc(wks_2,xpos1(i),0.865,mkres) gsn_text_ndc (wks_2,labels1(i),xpos2(i),0.865,txres) end do ;scatter_1 = gsn_csm_xy(wks_2,eof_ts(0,0:19),eof_ts(1,0:19),phase_res) scatter_1 = gsn_csm_xy(wks_2,eof_ts(0,0),eof_ts(1,0),phase_res) dum = new(4, graphic) ypts = (/1,1,-1,-1,1/) xpts = (/1,-1,-1,1,1/) do i = 0,3 dum(i) = gsn_add_polyline(wks_2,scatter_1,xpts(i:i+1),ypts(i:i+1),res_pyline) end do phase_res@xyMarkerColor = "green" ;scatter_2 = gsn_csm_xy(wks_2,eof_ts(0,20:39),eof_ts(1,20:39),phase_res) phase_res@xyMarkerColor = "blue" ;scatter_3 = gsn_csm_xy(wks_2,eof_ts(0,40:89),eof_ts(1,40:89),phase_res) ;overlay(scatter_1,scatter_2) ;overlay(scatter_1,scatter_3) dum_2 = new(69, graphic) do i = 1,69 if (i .le. 19) then mkres@gsMarkerColor = "green" end if if (i .gt. 19) then mkres@gsMarkerColor = "blue" end if dum_2(i-1) = gsn_add_polymarker(wks_2,scatter_1,eof_ts(0,i),eof_ts(1,i),mkres) end do mkres@gsMarkerIndex = 16 mkres@gsMarkerColor = "green" mkres@gsMarkerSizeF = 0.018 mkres@gsMarkerThicknessF = 10 mean_x = avg(eof_ts(0,0:19)) mean_y = avg(eof_ts(1,0:19)) gsn_polymarker(wks_2,scatter_1,mean_x,mean_y,mkres) mkres@gsMarkerColor = "blue" mean_x = avg(eof_ts(0,20:69)) mean_y = avg(eof_ts(1,20:69)) gsn_polymarker(wks_2,scatter_1,mean_x,mean_y,mkres) draw(scatter_1) frame(wks_2) ;;;;;;;;;;;;;;;;;;;;;;;;;;; Calculate the Cluster Means and Plot Them ;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;calculate_and_plot_cluster_forecasts(eof_ts,var,var_mean,lat_0,lon_0,ensemble0,"500 hPa Heights","HGT_500/cluster_hgt_500_forecasts.png") ;;;;;;;;;;;;;;;;;;;;;;; Calculate Cluster Images for Other Variables ;;;;;;;;;;;;;;;;;;;; max_temp_cluster_forecasts(eof_ts) min_temp_cluster_forecasts(eof_ts) qpf_cluster_forecasts(eof_ts) ;;;;;;;;;;;;;;;;;;;;;;; Calculate Cluster Evolution ;;;;;;;;;;;;;;;;;;;; ;calculate_and_plot_cluster_evolution(eof_ts) end