;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Author : Jeremy Berman; modifies by Kevin Biernat ; Created : 01/10/2017 ; Update : May 2019 ; Objective : PLOT COMPOSITE DIFFERENCES ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 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/contrib/ut_string.ncl" begin ;****************************************************************** ; set your mostAccurate/leastAccurate members ;****************************************************************** inds_mostAccurate = (/2,9,15,23,25,29,32,34,44,46/) inds_leastAccurate = (/4,8,19,21,27,30,31,41,45,48/) ;****************************************************************** ; Read in data ;****************************************************************** fileName_pf = "/keyserlab_rit/kbiernat/Crawford_cyclone_tracking/June_2018/composite_differences/AC1_108h_leadTime_init_2018053012/data/g_pf.grb2" inFile_pf = addfile( fileName_pf, "r" ) varin_pf = inFile_pf->gh_P1_L100_GLL0(:,:,{30000},:,:) fileName_cf = "/keyserlab_rit/kbiernat/Crawford_cyclone_tracking/June_2018/composite_differences/AC1_108h_leadTime_init_2018053012/data/g_cf.grb2" inFile_cf = addfile( fileName_cf, "r" ) varin_cf = inFile_cf->gh_P1_L100_GLL0(:,{30000},:,:) fileName_bounds = "/keyserlab_rit/kbiernat/Crawford_cyclone_tracking/June_2018/composite_differences/AC1_108h_leadTime_init_2018053012/bootstrap_bounds/g300_bootstrap_bounds.nc" ;-Output filename inFile_bounds = addfile( fileName_bounds, "r" ) mean_diff_lowerBound = inFile_bounds->mean_diff_lowerBound(:,:,:) mean_diff_upperBound = inFile_bounds->mean_diff_upperBound(:,:,:) forecast_hour = mean_diff_lowerBound&forecast_hour time = inFile_bounds->time ;********************************************************************************************** ; First combined all ensembles (control + perturbed) into one array and calculate ensemble mean ;********************************************************************************************** dimVar = dimsizes(varin_pf) nTimes = dimVar(1) nLat = dimVar(2) nLon = dimVar(3) var_all = new((/51,nTimes,nLat,nLon/),"float") var_all(0,:,:,:) = varin_cf var_all(1:50,:,:,:) = varin_pf var_all!0 = "ensemble0" var_all!1 = "forecast_hour" var_all!2 = "lat" var_all!3 = "lon" var_all&ensemble0 = ispan(0,50,1) var_all&forecast_hour = varin_pf&forecast_time0 var_all&lat = varin_pf&lat_0 var_all&lon = varin_pf&lon_0 var_avg = dim_avg_n_Wrap(var_all,0) ; ensemble mean var_avg = var_avg*0.1 ;************************************************************************************************************************ ; Subset data into composite groups, calculate statistical significance between groups ;************************************************************************************************************************ var_mostAccurate_4D = var_all(inds_mostAccurate,:,:,:) var_leastAccurate_4D = var_all(inds_leastAccurate,:,:,:) ; compute avg of mostAccurate and leastAccurate member groups var_avg_mostAccurate_3D = dim_avg_n_Wrap(var_mostAccurate_4D,0) var_avg_leastAccurate_3D = dim_avg_n_Wrap(var_leastAccurate_4D,0) ;************************************************************************************************************************ ; Calculate normalized composite differences ;************************************************************************************************************************ ; compute composite difference between groups var_diff_mostAccurate_leastAccurate_3D = var_avg_mostAccurate_3D - var_avg_leastAccurate_3D ; copy_VarMeta(dim_avg_n_Wrap(var_mostAccurate_4D,0),var_diff_mostAccurate_leastAccurate_3D) ; to copy coordinate variables to the composite difference variable var_std = var_diff_mostAccurate_leastAccurate_3D(:,:,:) ; to hold normalized composite difference ; below normalizes the field by the standard deviation; var_std_allmem = dim_stddev_n_Wrap(var_all, 0) ; ensemble standard deviation var_std = var_std/var_std_allmem ; divide by total ensemble standard deviation ;statSig = where(((var_diff_mostAccurate_leastAccurate_3D.lt.mean_diff_lowerBound).or.(var_diff_mostAccurate_leastAccurate_3D.gt.mean_diff_upperBound)),1.0,var_diff_mostAccurate_leastAccurate_3D@_FillValue) statSig = where(((var_diff_mostAccurate_leastAccurate_3D.lt.mean_diff_lowerBound).or.(var_diff_mostAccurate_leastAccurate_3D.gt.mean_diff_upperBound)),1.0,0.0) copy_VarMeta(var_diff_mostAccurate_leastAccurate_3D,statSig) ;print(statSig(0,:,:)) ;*********************************************************************************************************** ; Specify Attributes for Map Settings ;*********************************************************************************************************** res = True ; Turn on resources for map res@gsnDraw = False ; Do not draw plot until specified otherwise res@gsnFrame = False ; Do not advance frame until specified otherwise ;; Russia domain res@mpProjection = "Stereographic" res@mpCenterLatF = 90 res@mpCenterLonF = 60 res@mpGridPolarLonSpacingF = 0 res@mpLimitMode = "Corners" res@mpLeftCornerLatF = 20 res@mpRightCornerLatF = 59 res@mpLeftCornerLonF = 22 res@mpRightCornerLonF = 173 res@mpGridAndLimbOn = True res@mpGridLatSpacingF = 20 res@mpGridLonSpacingF = 20 ;res@mpCenterLonF = -40 res@mpGridLineColor = (/0.25,0.25,0.25/) res@mpGridLineThicknessF = 1.0 res@mpGridLineDashPattern = 10. res@mpGridAndLimbDrawOrder = "Draw" res@mpFillOn = True ; Turn off map fill ;res@mpOutlineDrawOrder = "Draw" res@mpFillDrawOrder = "PreDraw" res@mpOutlineOn = True ; Turn the map outline on res@mpOutlineBoundarySets = "National" ; draw national boundaries res@mpOutlineSpecifiers = (/"Canada", "United States"/) + " : States" ; draw US States, Canadian provinces res@mpOutlineDrawOrder = "Draw" res@mpGeophysicalLineThicknessF = 2.0 ; Change Thickness of the Geophysical Boundary Lines res@mpNationalLineThicknessF = 1.5 ; Change Thickness of National Boundary Lines res@mpUSStateLineThicknessF = 1.5 ; Change Thickness of US State and Canadian Province Boundaries res@mpGeophysicalLineColor = (/0.22,0.22,0.22/) res@mpNationalLineColor = res@mpGeophysicalLineColor res@mpUSStateLineColor = res@mpGeophysicalLineColor res@mpDataSetName = "Earth..4" ; Basically, this dataset draws more accurate coastlines/boundaries ; as well as state/provinces outlines. res@mpDataBaseVersion = 1 ; medium resolution res@mpDataResolution = 2 ; medium data resolution ;res@gsnMaximize = True ; make as large as possible res@gsnPaperOrientation = "landscape" ; Change orientation of display to landscape ;*********************************************************************************************************** ; Attributes for 500mb Height Contours ;*********************************************************************************************************** gres = True ; Resources for MSLP desired gres@gsnDraw = False ; Do not draw plot until specified otherwise gres@gsnFrame = False ; Do not advance frame until specified otherwise gres@gsnLeftString = "" ; Make blank the string above the plot to the left (usually a string representing variable name) gres@gsnRightString = "" ; Make blank the string above the plot to the right (usually a string representing variable units) gres@gsnAddCyclic = True gres@cnLinesOn = True ; Dispalay contour lines gres@cnLevelSelectionMode = "ExplicitLevels" ; use explicit levels gres@cnLevels = ispan(780,1200,12) ; set the contour levels gres@cnLineColor = 1 ; Color lines the foreground color (black) gres@cnFillOn = False ; Turn off contour fill. gres@cnLineLabelsOn = True ; Turn on line labels. gres@cnInfoLabelOn = False ; Turn off info label. gres@cnLineThicknessF = 5.0 gres@cnLabelMasking = True ; Contour lines are masked where contour labels appear gres@cnLineLabelInterval = 1 ; set label interval to every contour gres@cnLineLabelBackgroundColor = -1 ; set background box of label to transparent gres@cnMonoLineLabelFontColor = True ; Allow font colors to be displayed on all labels gres@cnLineLabelFontColor = gres@cnLineColor ; set color of labels on contours (can only use color table) gres@cnLineLabelFontThicknessF = 0.0 ; Line Label Font thickness (multiplier of unit thickness of 1.0) gres@cnLineLabelFontHeightF = 0.014 ; Change Font Size of Contour Labels gres@cnHighLabelPerimSpaceF = 0.00 gres@cnLineLabelDensityF = 1.0 gres@cnLineDrawOrder = "Draw" gres@cnLabelDrawOrder = "Draw" gres@cnLineLabelAngleF = 0.0 gres@cnLineLabelPlacementMode = "Computed" ;************************************************ ; Resources for composite differences ;************************************************ fres =True fres@gsnAddCyclic = True fres@gsnDraw = False fres@gsnFrame = False fres@gsnRightString ="" fres@gsnLeftString = "" fres@cnFillOn = True fres@cnFillDrawOrder = "PreDraw" fres@cnLinesOn = False fres@cnLineThicknessF = .5 fres@cnLineLabelsOn = False fres@gsnSpreadColors = False fres@cnLevelSelectionMode = "ExplicitLevels" ;fres@cnLevels = (/-2.0,-1.75,-1.5,-1.25,-1.0,-0.75,-0.5,-0.25,0,0.25,0.5,0.75,1.0,1.25,1.5,1.75,2.0/) ;fres@cnFillColors = (/10,20,30,40,50,60,70,85,-1,-1,170,185,195,205,215,225,235,245/) fres@cnLevels = (/-2.0,-1.6,-1.2,-0.8,-0.4,0,0.4,0.8,1.2,1.6,2.0/) fres@cnFillColors = (/12,28,43,58,85,-1,-1,170,190,208,225,250/) fres@cnSmoothingOn = True fres@cnSmoothingTensionF = -2. fres@cnSmoothingDistanceF = 0.001 fres@lbLabelBarOn = True fres@cnInfoLabelOn = False fres@lbOrientation = "Horizontal" fres@pmLabelBarOrthogonalPosF = 0.00 fres@lbLabelStride = 1 fres@pmLabelBarHeightF = 0.06 fres@lbLabelFontHeightF = .010 fres@lbBoxLinesOn = True ;*********************************************************************************************************** ; Attributes for 95% confidence level contour ;*********************************************************************************************************** scres = True scres@gsnDraw = False ; Do not draw plot until specified otherwise scres@gsnFrame = False ; Do not advance frame until specified otherwise scres@gsnLeftString = "" scres@gsnRightString = "" scres@gsnAddCyclic = True scres@cnLinesOn = True ; Dispalay contour lines scres@cnLevelSelectionMode = "ExplicitLevels" scres@cnLevels = 0.9 ; set the minimum contour level scres@cnLineColor = (/1.00,1.00,1.00/) ; scres@cnFillOn = False ; Turn off contour fill. scres@cnLineLabelsOn = False ; Turn on line labels. scres@cnInfoLabelOn = False ; Turn off info label. scres@cnLineThicknessF = 4.0 scres@cnLabelMasking = True ; Contour lines are masked where contour labels appear scres@cnLineLabelDensityF = 1.0 scres@cnLineLabelAngleF = 0.0 ;scres@cnLineLabelPlacementMode = 2 scres@cnLineLabelBackgroundColor = -1 ; set background box of label to transparent scres@cnMonoLineLabelFontColor = True ; Allow font colors to be displayed on all labels ;scres@cnLineLabelFont = "Helvetica" scres@cnLineLabelInterval = 1 scres@cnLineLabelFontColor = scres@cnLineColor ; set color of labels on contours (can only use color table) scres@cnLineLabelFontThicknessF = 1.5 ; Line Label Font thickness (multiplier of unit thickness of 1.0) scres@cnLineLabelFontHeightF = 0.016 ; Change Font Size of Contour Labels scres@cnLineLabelPerimSpaceF = 0.00 ; Essentially, limit the margins surrounding each label bo scres@cnLineLabelPlacementMode = "Computed" scres@cnLineDrawOrder = "Draw" scres@cnLabelDrawOrder = "Draw" scres@cnSmoothingOn = True scres@cnSmoothingTensionF = -2. scres@cnSmoothingDistanceF = 0.001 ;*********************************************************************************************************** ; Attributes for shading above 95% confidence level ;*********************************************************************************************************** sfres =True sfres@gsnAddCyclic = True sfres@gsnDraw = False sfres@gsnFrame = False sfres@gsnRightString ="" sfres@gsnLeftString = "" sfres@cnFillOn = True sfres@cnFillPattern = (/17/) sfres@cnFillDrawOrder = "Draw" sfres@cnLinesOn = False sfres@cnLineThicknessF = .5 sfres@cnLineLabelsOn = False sfres@gsnSpreadColors = False sfres@cnLevelSelectionMode = "ExplicitLevels" sfres@cnLevels = (/0.9/) sfres@cnFillColors = (/-1,0/) sfres@cnFillDotSizeF = 0.0021 sfres@cnSmoothingOn = True sfres@cnSmoothingTensionF = -2. sfres@cnSmoothingDistanceF = 0.002 sfres@lbLabelBarOn = False sfres@cnInfoLabelOn = False sfres@lbOrientation = "Horizontal" sfres@pmLabelBarOrthogonalPosF = 0.00 sfres@lbLabelStride = 1 sfres@pmLabelBarHeightF = 0.05 sfres@lbLabelFontHeightF = .010 sfres@lbBoxLinesOn = True ;****************************************************************** ; Plotting ;****************************************************************** do i=0, nTimes-1 wks_type = "png" wks_type@wkWidth = 1536 wks_type@wkHeight = 1536 imageDirectory = "/keyserlab_rit/kbiernat/Crawford_cyclone_tracking/June_2018/composite_differences/AC1_108h_leadTime_init_2018053012/images/g300" imageName = "compd_g300_F"+forecast_hour(i)+"_valid_"+ut_string(time(i),"%Y%N%D_%H%M") wks = gsn_open_wks(wks_type, imageName) ; open a ps file with a unique name according to the time/date gsn_define_colormap(wks,"BlueWhiteOrangeRed") res@gsnLeftString = ut_string(time(i),"%H%M UTC %d %c %Y") res@gsnRightString = "" res@gsnLeftStringFontHeightF = 0.01 res@gsnRightStringFontHeightF = 0.0 plot_map = gsn_csm_map(wks,res) plot_compd = gsn_csm_contour(wks,var_std(i,:,:),fres) plot_sigf = gsn_csm_contour(wks,statSig(i,:,:),sfres) plot_sigc = gsn_csm_contour(wks,statSig(i,:,:),scres) plot_mean = gsn_csm_contour(wks,var_avg(i,:,:),gres) overlay(plot_map,plot_compd) overlay(plot_map,plot_sigf) overlay(plot_map,plot_sigc) overlay(plot_map,plot_mean) res@gsnDraw = False ; Do not draw plot until specified otherwise res@gsnFrame = False ; Do not advance frame until specified otherwise maximize_output(wks,res) draw(plot_map) ; Finally, draw the Plot frame( wks ) ; Finally, advance the frame system("convert -trim +repage "+imageName+".png "+ \ imageDirectory+"/"+imageName+".png") system("rm *.png") end do end