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 procedure adds metadata (lat, lon, ensemble) to a 2 dimensional variable. undef("createVarMeta_2d") procedure createVarMeta_2d(x:numeric, lat[*]:numeric, lon[*]:numeric) ; add meta data begin x@_FillValue = 1e20 dimx = dimsizes(x) rankx = dimsizes(dimx) if (rankx.eq.2) then x!0 = "lat_0" x!1 = "lon_0" x&lat_0 = lat x&lon_0 = lon else print("createVarMeta: error: rankx="+rankx) exit end if end ; This function finds the corresponding name for the input numeric month e.g., 04 = Apr undef("monthName") function monthName(numeric_month) begin month_name_text = "Jan" if(numeric_month .eq. 1) then month_name_text = "Jan" end if if(numeric_month .eq. 2) then month_name_text = "Feb" end if if(numeric_month .eq. 3) then month_name_text = "Mar" end if if(numeric_month .eq. 4) then month_name_text = "Apr" end if if(numeric_month .eq. 5) then month_name_text = "May" end if if(numeric_month .eq. 6) then month_name_text = "Jun" end if if(numeric_month .eq. 7) then month_name_text = "Jul" end if if(numeric_month .eq. 8) then month_name_text = "Aug" end if if(numeric_month .eq. 9) then month_name_text = "Sep" end if if(numeric_month .eq. 10) then month_name_text = "Oct" end if if(numeric_month .eq. 11) then month_name_text = "Nov" end if if(numeric_month .eq. 12) then month_name_text = "Dec" end if return(month_name_text) 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 = tointeger(sum(cluster_mem(0:19))) cmce_percentage = round((tofloat(cmce_membership)/20.0)*100,3) gefs_membership = tointeger(sum(cluster_mem(20:39))) gefs_percentage = round((tofloat(gefs_membership)/20.0)*100,3) ens_membership = tointeger(sum(cluster_mem(40:89))) ens_percentage = round((tofloat(ens_membership)/50.0)*100,3) total_membership = tointeger(cmce_membership+gefs_membership+ens_membership) total_percentage = round((tofloat(total_membership)/90.0)*100,3) return([/cmce_membership,cmce_percentage,gefs_membership,gefs_percentage,ens_membership,ens_percentage,total_membership,total_percentage/]) end ; This fuction calculates and plots the cluster forecasts undef("calculate_and_plot_cluster_forecasts") procedure calculate_and_plot_cluster_forecasts(eof_ts, clusters, var, var_mean, lat_0, lon_0, ensemble0, var_name, cluster_image_name, init_year_str, init_month_txt, init_day_str) begin cluster_1_members = where(clusters@id .eq. 1,1.0,eof_ts@_FillValue) cluster_1_members_conform = conform(var, cluster_1_members, (/2/)) cluster_1 = where(cluster_1_members_conform .eq. 1.0, var, var@_FillValue) createVarMeta(cluster_1, lat_0, lon_0, ensemble0) cluster_1_mean = dim_avg_Wrap(cluster_1(lat_0|:,lon_0|:,ensemble0|:)) cluster_1_diff = cluster_1_mean - var_mean copy_VarMeta(cluster_1_mean, cluster_1_diff) cluster_2_members = where(clusters@id .eq. 2,1.0,eof_ts@_FillValue) cluster_2_members_conform = conform(var, cluster_2_members, (/2/)) cluster_2 = where(cluster_2_members_conform .eq. 1.0, var, var@_FillValue) createVarMeta(cluster_2, lat_0, lon_0, ensemble0) cluster_2_mean = dim_avg_Wrap(cluster_2(lat_0|:,lon_0|:,ensemble0|:)) cluster_2_diff = cluster_2_mean - var_mean copy_VarMeta(cluster_2_mean, cluster_2_diff) cluster_3_members = where(clusters@id .eq. 3,1.0,eof_ts@_FillValue) cluster_3_members_conform = conform(var, cluster_3_members, (/2/)) cluster_3 = where(cluster_3_members_conform .eq. 1.0, var, var@_FillValue) createVarMeta(cluster_3, lat_0, lon_0, ensemble0) cluster_3_mean = dim_avg_Wrap(cluster_3(lat_0|:,lon_0|:,ensemble0|:)) cluster_3_diff = cluster_3_mean - var_mean copy_VarMeta(cluster_3_mean, cluster_3_diff) cluster_4_members = where(clusters@id .eq. 4,1.0,eof_ts@_FillValue) cluster_4_members_conform = conform(var, cluster_4_members, (/2/)) cluster_4 = where(cluster_4_members_conform .eq. 1.0, var, var@_FillValue) createVarMeta(cluster_4, lat_0, lon_0, ensemble0) cluster_4_mean = dim_avg_Wrap(cluster_4(lat_0|:,lon_0|:,ensemble0|:)) cluster_4_diff = cluster_4_mean - var_mean copy_VarMeta(cluster_4_mean, cluster_4_diff) total_mean = dim_avg_Wrap(var(lat_0|:,lon_0|:,ensemble0|:)) total_mean_diff = total_mean - var_mean copy_VarMeta(total_mean, total_mean_diff) total_mean_diff(0,0) = total_mean_diff(0,0)+0.001 imagedir = "/cf_hmt/hmt/lamberson/day_8_10/fuzzy_clustering/images/cluster_forecast/" cluster_wks_name = imagedir+cluster_image_name wks_type = "png" wks_type@wkWidth = 1800 wks_type@wkHeight = 1800 cluster_wks = gsn_open_wks(wks_type,cluster_wks_name) if (var_name .eq. "qpf") then color_table_name="precip4_diff_19lev" fill_colors = (/2,3,5,6,7,8,9,10,0,13,14,15,16,17,18,19,21/) else color_table_name="amwg_blueyellowred" fill_colors = (/2,3,4,5,6,7,8,9,0,10,11,12,13,14,15,16,17/) end if gsn_define_colormap(cluster_wks,color_table_name) cluster_diff_res = True cluster_diff_res@gsnDraw = False cluster_diff_res@gsnFrame = False cluster_diff_res@mpProjection = "CylindricalEquidistant" ;cluster_diff_res@gsnMaskLambertConformal = True cluster_diff_res@mpDataBaseVersion = "MediumRes" cluster_diff_res@mpLimitMode = "LatLon" cluster_diff_res@mpMinLatF = 25. cluster_diff_res@mpMaxLatF = 50. cluster_diff_res@mpMaxLonF = 285. cluster_diff_res@mpMinLonF = 245. 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 = (/-240.,-210.,-180.,-150.,-120.,-90.,-60.,-30.,30.,60.,90.,120.,150.,180.,210.,240./) cluster_diff_res@cnFillOn = True cluster_diff_res@cnFillColors = fill_colors 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.017 cluster_diff_res@gsnRightStringFontHeightF = 0.017 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 = 4 cluster_mean_res@gsnContourNegLineDashPattern = 1 if(var_name .eq. "500 hPa Heights") then cluster_diff_res@mpMinLatF = 25. cluster_diff_res@mpMaxLatF = 80. cluster_diff_res@mpMaxLonF = 300. cluster_diff_res@mpMinLonF = 225. 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(cluster_1_members) cluster_diff_res@gsnLeftString = "Cluster 1" cluster_diff_res@gsnRightString = "C: "+tostring(membership[0])+"="+tostring(membership[1])+"% G: "+tostring(membership[2])+"="+tostring(membership[3])+"% E: "+tostring(membership[4])+"="+tostring(membership[5])+"% T: "+tostring(membership[6])+"="+tostring(membership[7])+"%" cluster_diff_res@gsnRightStringFontHeightF=0.017 cluster_diff_res@gsnLeftStringFontHeightF=0.017 cluster_diff_res@gsnCenterString = "" cluster_1_diff_plot = gsn_csm_contour_map(cluster_wks,cluster_1_diff,cluster_diff_res) cluster_1_mean_plot = gsn_csm_contour(cluster_wks,cluster_1_mean,cluster_mean_res) overlay(cluster_1_diff_plot,cluster_1_mean_plot) membership = clusterMembershipCount(cluster_2_members) cluster_diff_res@gsnLeftString = "Cluster 2" cluster_diff_res@gsnRightString = "C: "+tostring(membership[0])+"="+tostring(membership[1])+"% G: "+tostring(membership[2])+"="+tostring(membership[3])+"% E: "+tostring(membership[4])+"="+tostring(membership[5])+"% T: "+tostring(membership[6])+"="+tostring(membership[7])+"%" cluster_diff_res@gsnCenterString = "" cluster_2_diff_plot = gsn_csm_contour_map(cluster_wks,cluster_2_diff,cluster_diff_res) cluster_2_mean_plot = gsn_csm_contour(cluster_wks,cluster_2_mean,cluster_mean_res) overlay(cluster_2_diff_plot,cluster_2_mean_plot) membership = clusterMembershipCount(cluster_3_members) cluster_diff_res@gsnLeftString = "Cluster 3" cluster_diff_res@gsnRightString = "C: "+tostring(membership[0])+"="+tostring(membership[1])+"% G: "+tostring(membership[2])+"="+tostring(membership[3])+"% E: "+tostring(membership[4])+"="+tostring(membership[5])+"% T: "+tostring(membership[6])+"="+tostring(membership[7])+"%" cluster_diff_res@gsnCenterString = "" cluster_3_diff_plot = gsn_csm_contour_map(cluster_wks,cluster_3_diff,cluster_diff_res) cluster_3_mean_plot = gsn_csm_contour(cluster_wks,cluster_3_mean,cluster_mean_res) overlay(cluster_3_diff_plot,cluster_3_mean_plot) membership = clusterMembershipCount(cluster_4_members) cluster_diff_res@gsnLeftString = "Cluster 4" cluster_diff_res@gsnRightString = "C: "+tostring(membership[0])+"="+tostring(membership[1])+"% G: "+tostring(membership[2])+"="+tostring(membership[3])+"% E: "+tostring(membership[4])+"="+tostring(membership[5])+"% T: "+tostring(membership[6])+"="+tostring(membership[7])+"%" cluster_diff_res@gsnCenterString = "" cluster_4_diff_plot = gsn_csm_contour_map(cluster_wks,cluster_4_diff,cluster_diff_res) cluster_4_mean_plot = gsn_csm_contour(cluster_wks,cluster_4_mean,cluster_mean_res) overlay(cluster_4_diff_plot,cluster_4_mean_plot) cluster_diff_res@gsnLeftString = "Total Mean" cluster_diff_res@gsnRightString = "C: 20=100% G: 20=100% E: 50=100% T: 90=100%" cluster_diff_res@gsnCenterString = "" total_mean_diff_plot = gsn_csm_contour_map(cluster_wks,total_mean_diff,cluster_diff_res) total_mean_plot = gsn_csm_contour(cluster_wks,total_mean,cluster_mean_res) overlay(total_mean_diff_plot,total_mean_plot) cluster_hgt_plots = new(5,graphic) cluster_hgt_plots(0) = cluster_1_diff_plot cluster_hgt_plots(1) = cluster_2_diff_plot cluster_hgt_plots(2) = cluster_3_diff_plot cluster_hgt_plots(3) = cluster_4_diff_plot cluster_hgt_plots(4) = total_mean_diff_plot panel_res = True panel_res@gsnMaximize = True panel_res@txString = "Init: 0000 UTC "+init_month_txt+" "+init_day_str+" "+init_year_str panel_res@txFont = 22 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,(/3,2/),panel_res) delete(cluster_wks) delete(cluster_hgt_plots) cmd = "convert -trim -geometry 1800x1800 +repage -border 8 -bordercolor white -background white -flatten "+cluster_wks_name+" "+cluster_wks_name system(cmd) end ; This procedure calculates and plots the cluster forecasts of maximum temperature undef("max_temp_cluster_forecasts") procedure max_temp_cluster_forecasts(eof_ts, clusters, init_year_str, init_month_txt, init_day_str) 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,89,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(cmce_max_temps,gefs_max_temps,0) max_temp_combined_2 = array_append_record(max_temp_combined_1,ens_max_temps,0) max_temp = max_temp_combined_2(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_k_means_central.png" max_temp_3_day_sum = max_temp else if(i .eq. 1) then cluster_image_name = "TMAX/day_9_cluster_maximum_temperatures_k_means_central.png" max_temp_3_day_sum = max_temp_3_day_sum+max_temp else cluster_image_name = "TMAX/day_10_cluster_maximum_temperatures_k_means_central.png" max_temp_3_day_sum = max_temp_3_day_sum+max_temp end if end if calculate_and_plot_cluster_forecasts(eof_ts,clusters,max_temp,max_temp_mean,lat_0,lon_0,ensemble0,"maximum temperature",cluster_image_name, init_year_str, init_month_txt, init_day_str) j = j+4 delete([/max_temp,cmce_max_temps,gefs_max_temps,ens_max_temps,max_temp_combined_1,max_temp_combined_2,max_temp_mean/]) end do max_temp_3_day_mean = max_temp_3_day_sum/3.0 createVarMeta(max_temp_3_day_mean, lat_0, lon_0, ensemble0) max_temp_3_day_total_mean = dim_avg_Wrap(max_temp_3_day_mean(lat_0|:,lon_0|:,ensemble0|:)) cluster_image_name = "TMAX/3_day_avg_cluster_maximum_temperatures_k_means_central.png" calculate_and_plot_cluster_forecasts(eof_ts,clusters,max_temp_3_day_mean,max_temp_3_day_total_mean,lat_0,lon_0,ensemble0,"maximum temperature",cluster_image_name, init_year_str, init_month_txt, init_day_str) end undef("min_temp_cluster_forecasts") procedure min_temp_cluster_forecasts(eof_ts, clusters, init_year_str, init_month_txt, init_day_str) 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,89,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(cmce_min_temps,gefs_min_temps,0) min_temp_combined_2 = array_append_record(min_temp_combined_1,ens_min_temps,0) min_temp = min_temp_combined_2(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_k_means_central.png" min_temp_3_day_sum = min_temp else if(i .eq. 1) then cluster_image_name = "TMIN/day_9_cluster_minimum_temperatures_k_means_central.png" min_temp_3_day_sum = min_temp_3_day_sum+min_temp else cluster_image_name = "TMIN/day_10_cluster_minimum_temperatures_k_means_central.png" min_temp_3_day_sum = min_temp_3_day_sum+min_temp end if end if calculate_and_plot_cluster_forecasts(eof_ts,clusters,min_temp,min_temp_mean,lat_0,lon_0,ensemble0,"minimum temperature",cluster_image_name, init_year_str, init_month_txt, init_day_str) j = j+4 delete([/min_temp,cmce_min_temps,gefs_min_temps,ens_min_temps,min_temp_combined_1,min_temp_combined_2,min_temp_mean/]) end do min_temp_3_day_mean = min_temp_3_day_sum/3.0 createVarMeta(min_temp_3_day_mean, lat_0, lon_0, ensemble0) min_temp_3_day_total_mean = dim_avg_Wrap(min_temp_3_day_mean(lat_0|:,lon_0|:,ensemble0|:)) cluster_image_name = "TMIN/3_day_avg_cluster_minimum_temperatures_k_means_central.png" calculate_and_plot_cluster_forecasts(eof_ts,clusters,min_temp_3_day_mean,min_temp_3_day_total_mean,lat_0,lon_0,ensemble0,"minimum temperature",cluster_image_name, init_year_str, init_month_txt, init_day_str) end undef("qpf_cluster_forecasts") procedure qpf_cluster_forecasts(eof_ts, clusters, init_year_str, init_month_txt, init_day_str) 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,89,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(cmce_qpf,gefs_qpf,0) qpf_combined_2 = array_append_record(qpf_combined_1,ens_qpf,0) qpf = qpf_combined_2(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_k_means_central.png" calculate_and_plot_cluster_forecasts(eof_ts,clusters,qpf,qpf_mean,lat_0,lon_0,ensemble0,"qpf",cluster_image_name, init_year_str, init_month_txt, init_day_str) delete([/qpf,cmce_qpf,gefs_qpf,ens_qpf,qpf_combined_1,qpf_combined_2,qpf_mean/]) 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) 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(:,:,2,:,:) ens_var_4d_evolution = ens_var_4d_evolution_temp(ensemble0|:,ncl_join|:,g0_lat_2|:,g0_lon_3|:) 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) 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 = "+ EOF1 Cluster" cluster_diff_res@gsnRightString = "" cluster_diff_res@gsnCenterString = "" 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,"WhiteBlueGreenYellowRed") 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,"WhiteBlueGreenYellowRed") cluster_diff_res@gsnLeftString = "- 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,"WhiteBlueGreenYellowRed") cluster_diff_res@gsnLeftString = "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,"WhiteBlueGreenYellowRed") cluster_diff_res@gsnLeftString = "+ 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,"WhiteBlueGreenYellowRed") cluster_diff_res@gsnLeftString = "- 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 ;;;;;;;;;;;;; Format the initialization date ;;;;;;;;;;;;;;;;;; init_year_str = tostring(init_year) init_month_txt = monthName(init_month) if(init_month .lt. 10) then init_month_str = "0"+tostring(init_month) else init_month_str = tostring(init_month) end if if(init_day .lt. 10) then init_day_str = "0"+tostring(init_day) else init_day_str = tostring(init_day) end if ;;;;;;;;;;;;;; Format the forecast date ;;;;;;;;;;;;;;;;;;;;;; fcst_year_str = tostring(fcst_year) fcst_month_txt = monthName(fcst_month) if(fcst_month .lt. 10) then fcst_month_str = "0"+tostring(fcst_month) else fcst_month_str = tostring(fcst_month) end if if(fcst_day .lt. 10) then fcst_day_str = "0"+tostring(fcst_day) else fcst_day_str = tostring(fcst_day) end if 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 = dim_avg_n_Wrap(ens_var_4d,(/0/)) 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(cmce_var,gefs_var,0) dummy_combined_var_2 = array_append_record(dummy_combined_var_1,ens_var,0) var = dummy_combined_var_2(dim1|:,dim2|:,dim0|:) lat_0 = gefs_f[0]->lat_0 lon_0 = gefs_f[0]->lon_0 ensemble0 = ispan(0,89,1) createVarMeta(var, lat_0, lon_0, ensemble0) var_mean = dim_avg_Wrap(var(lat_0|:,lon_0|:,ensemble0|:)) var_mean_conform = conform(var, var_mean,(/0,1/)) 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 = 20.0 max_lon = 300.0 min_lon = 225.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) opt = True opt@iseed = 2 clusters = kmeans_as136(eof_ts,4,opt) wks_name=imagedir+"EOFs_k_means_central.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@txString = "Init: 0000 UTC "+init_month_txt+" "+init_day_str+" "+init_year_str panel_res@txFont = 22 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,(/2/),panel_res) wks_name_2=imagedir+"cluster_phase_space_k_means_central.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 = "Init: 0000 UTC "+init_month_txt+" "+init_day_str+" "+init_year_str phase_res@tiMainFont = 22 phase_res@tiMainString = "Scatter Plot" phase_res@xyMarkLineMode = "Markers" phase_res@xyMarkers = 5 phase_res@xyMarkerColor = "red" phase_res@xyMarkerSizeF = 0.001 phase_res@xyMarkerThicknessF = 0.001 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 scatter_1 = gsn_csm_xy(wks_2,5,5,phase_res) mkres = True ; Marker resources mkres@gsMarkerSizeF = 0.015 mkres@gsMarkerThicknessF = 6.0 dum_2 = new(90, graphic) do i = 0,89 if (i .le. 19) then mkres@gsMarkerColor = "red" end if if (i .gt. 19 .and. i .lt. 40) then mkres@gsMarkerColor = "green" end if if (i .gt. 39) then mkres@gsMarkerColor = "blue" end if if(clusters@id(i) .eq. 1) then mkres @gsMarkerIndex = 5 end if if(clusters@id(i) .eq. 2) then mkres @gsMarkerIndex = 15 end if if(clusters@id(i) .eq. 3) then mkres @gsMarkerIndex = 9 end if if(clusters@id(i) .eq. 4) then mkres @gsMarkerIndex = 7 end if if(clusters@id(i) .eq. 5) then mkres @gsMarkerIndex = 6 end if dum_2(i) = gsn_add_polymarker(wks_2,scatter_1,eof_ts(0,i),eof_ts(1,i),mkres) end do mkres@gsMarkerIndex = 16 mkres@gsMarkerColor = "red" 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 = "green" mean_x = avg(eof_ts(0,20:39)) mean_y = avg(eof_ts(1,20:39)) gsn_polymarker(wks_2,scatter_1,mean_x,mean_y,mkres) mkres@gsMarkerColor = "blue" mean_x = avg(eof_ts(0,40:89)) mean_y = avg(eof_ts(1,40:89)) gsn_polymarker(wks_2,scatter_1,mean_x,mean_y,mkres) draw(scatter_1) ;---------------------------------------------------------------------- ; Draw some individual labelbars. ;---------------------------------------------------------------------- labels1 = (/"CMCE","GEFS","ECMWF"/) colors1 = (/"red","green","blue"/) lbres = True ; labelbar only resources lbres@vpWidthF = 0.1 ; labelbar width lbres@vpHeightF = 0.1 ; labelbar height lbres@lbBoxMajorExtentF = 0.15 ; puts space between color boxes lbres@lbMonoFillPattern = True ; Solid fill pattern lbres@lbLabelFontHeightF = 0.015 ; font height. default is small lbres@lbLabelJust = "CenterLeft" ; left justify labels lbres@lbPerimOn = False ; ; Each labelbar has just one label. This allows you to more ; easily control where the label goes. ; xpos = (/0.35, 0.50, 0.65/) do i=0,2 lbres@lbFillColors = colors1(i) lbres@lbLabelFontColor = colors1(i) gsn_labelbar_ndc(wks_2,1,labels1(i),xpos(i),0.93,lbres) end do ;---------------------------------------------------------------------- ; Draw some markers and text. ;---------------------------------------------------------------------- labels2 = (/"Cluster 1","Cluster 2","Cluster 3","Cluster 4"/) markers1 = (/ 5, 15, 9, 7/) xpos2 = (/ 0.21, 0.42, 0.63, 0.83/) xpos3 = (/ 0.23, 0.44, 0.65, 0.85/) mkres2 = True ; Marker resources txres = True ; Text resources txres@txFontHeightF = 0.015 txres@txJust = "CenterLeft" do i = 0,3 mkres2@gsMarkerThicknessF = 6.0 mkres2@gsMarkerSizeF = 0.015 mkres2@gsMarkerIndex = markers1(i) gsn_polymarker_ndc(wks_2, xpos2(i),0.15,mkres2) gsn_text_ndc (wks_2,labels2(i),xpos3(i),0.15,txres) end do frame(wks_2) ;;;;;;;;;;;;;;;;;;;;;;;;;;; Calculate the Cluster Means and Plot Them ;;;;;;;;;;;;;;;;;;;;;;;;;;;;; cfsr_dir = "/cf_hmt/hmt/lamberson/day_8_10/dprogdt/cfsr/" cfsr_file = cfsr_dir+"cfsrmean_"+fcst_month_str+"_"+fcst_day_str+"_0000Z.grb2" cfsr_f = addfile(cfsr_file, "r") cfsr_var = cfsr_f->HGT_P0_L100_GLL0({50000},:,:) ;cfsr_var_reduced = cfsr_var(::-1,::-1) calculate_and_plot_cluster_forecasts(eof_ts,clusters,var({min_lat:max_lat},{min_lon:max_lon},:),cfsr_var({min_lat:max_lat},{min_lon:max_lon}),lat_0({min_lat:max_lat}),lon_0({min_lon:max_lon}),ensemble0,"500 hPa Heights","HGT_500/cluster_hgt_500_forecasts_k_means_central.png",init_year_str,init_month_txt,init_day_str) ;;;;;;;;;;;;;;;;;;;;;;; Calculate Cluster Images for Other Variables ;;;;;;;;;;;;;;;;;;;; max_temp_cluster_forecasts(eof_ts,clusters,init_year_str,init_month_txt,init_day_str) min_temp_cluster_forecasts(eof_ts,clusters,init_year_str,init_month_txt,init_day_str) qpf_cluster_forecasts(eof_ts,clusters,init_year_str,init_month_txt,init_day_str) ;;;;;;;;;;;;;;;;;;;;;;; Calculate Cluster Evolution ;;;;;;;;;;;;;;;;;;;; ;calculate_and_plot_cluster_evolution(eof_ts) end