;*** Create maps for SSW outlook**** 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" begin ;*************USER DEFINED VARIABLES***************** ; define hail threshold size nsize = 1.0 ; define year through which composites created yyyy = 2017 ; define season of interest seas = "MAM" mo1 = 3 mo2 = 4 mo3 = 5 ; Define ONI thresholds for calculating anomaly maps range_values = new(2,"double") range_values(0) = 0.5 range_values(1) = 1.0 ; Define region over which count is performed *** latS = 25 latN = 50 lonW = -125 lonE = -65 ; Define file names and locations file_SST = "/export/cpc-lw-kharnos/kirstin.harnos/SSW/SST/oni.ascii.txt" file_hail = "/export/cpc-lw-kharnos/kirstin.harnos/SSW/HAIL/1955-2017_hail.csv" file_unique_hail = "/export/cpc-lw-kharnos/kirstin.harnos/SSW/HAIL/Unique_hail_days_1980-2017.nc" ; Define output image type and filename image_type = "png" out_fig_name = "outlook2019.MAM-hail-reports.MAM-ONI.composites" ;***********END USER DEFINE VARIABLES********************* ; Read ONI File ncol = 4 data = asciiread(file_SST,-1,"string") delim = " " nfields = str_fields_count(data(0),delim) season_sst = str_get_field(data,1,delim) year_sst = stringtofloat(str_get_field(data,2,delim)) nino_sst = stringtofloat(str_get_field(data,3,delim)) anom_sst = stringtofloat(str_get_field(data,4,delim)) ; define a time array and set as corrdinate variable time_sst = yyyymm_time(1950,yyyy,"integer") nino_sst!0 = "time" nino_sst&time = time_sst anom_sst!0 = "time" anom_sst&time = time_sst ; Filter data to 1980-yyyy and NDJ iyear = ind(year_sst .ge.1980 .and. year_sst.le.yyyy) year = year_sst(iyear) season = season_sst(iyear) nino = nino_sst(iyear) anom = anom_sst(iyear) delete(iyear) delete(year_sst) delete(season_sst) delete(nino_sst) delete(anom_sst) iseas = ind(season .eq. seas) year_sst = year(iseas) season_sst = season(iseas) nino_sst = nino(iseas) anom_sst1 = anom(iseas) anom_sst = dim_standardize_Wrap(anom_sst1,1) delete(year) delete(season) delete(nino) delete(anom) ;** Read Hail Counts ** filename = asciiread(file_hail,-1,"string") delim = "," ; Get desired fields year = tointeger(str_get_field(filename,2,delim)) month = tointeger(str_get_field(filename,3,delim)) lat = tofloat(str_get_field(filename,16,delim)) lon = tofloat(str_get_field(filename,17,delim)) size = tofloat(str_get_field(filename,11,delim)) delete(filename) ; Filter data to 'since 1980', MAM and CONUS iyear = ind(year.ge.1980) year_filter1 = year(iyear) month_filter1 = month(iyear) lat_filter1 = lat(iyear) lon_filter1 = lon(iyear) size_filter1 = size(iyear) delete(year) delete(month) delete(lat) delete(lon) delete(size) ifilter = ind(month_filter1.eq.mo1 .or. month_filter1.eq.mo2 .or. month_filter1.eq.mo3) year_mam = year_filter1(ifilter) month_mam = month_filter1(ifilter) size_mam = size_filter1(ifilter) lat_mam = lat_filter1(ifilter) lon_mam = lon_filter1(ifilter) ilon_lat = ind(lat_mam.le.latN .and. lat_mam.ge.latS .and. lon_mam.le.lonE .and. lon_mam.ge.lonW) year_CONUS = year_mam(ilon_lat) size_CONUS = size_mam(ilon_lat) lat_CONUS = lat_mam(ilon_lat) lon_CONUS = lon_mam(ilon_lat) delete(year_filter1) delete(size_filter1) delete(lat_filter1) delete(lon_filter1) delete(ifilter) ; Filter data to hail size isize = ind(size_CONUS.ge.nsize) year_reports = year_CONUS(isize) size_reports = size_CONUS(isize) lat_reports = lat_CONUS(isize) lon_reports = lon_CONUS(isize) delete(year_mam) delete(size_mam) delete(lat_mam) delete(lon_mam) delete(isize) delete(year_CONUS) delete(size_CONUS) delete(lat_CONUS) delete(lon_CONUS) delete(ilon_lat) ;** Read Hail Days data (1980-2014) created using MATLAB ** f = addfile(file_unique_hail,"r") ; Get desired fields year_Hdays = f->unique_year month_Hdays = f->unique_month lat_Hdays = f->unique_lat lon_Hdays = f->unique_lon size_Hdays = f->unique_size delete(f) ;*** Grid tornado/hail data **** ;--- Create desired grid dlat = 1 dlon = 1 nlat = toint((latN-latS)/dlat) mlon = toint((lonE-lonW)/dlon) lat = fspan(latS, latN, nlat) lon = fspan(lonW, lonE, mlon) lat@units = "degrees_north" lon@units = "degrees_east" nyears = (max(year_reports) - min(year_reports))+1 years = ispan(min(year_reports),max(year_reports),1) ;--- Bin data reports = new( (/nyears,nlat,mlon/), "float") reports!0 = "years" reports!1 = "lat" reports!2 = "lon" reports&years = years reports&lat = lat reports&lon = lon Hdays = new( (/nyears,nlat,mlon/), "float") Hdays!0 = "years" Hdays!1 = "lat" Hdays!2 = "lon" Hdays&years = years Hdays&lat = lat Hdays&lon = lon reports = 0 Hdays = 0 ; Grid Reports data do i = 0,dimsizes(years)-1 do n=0,dimsizes(size_reports)-1 if (year_reports(n) .eq. years(i)) then jl = toint((lat_reports(n)-latS)/dlat) il = toint((lon_reports(n)-lonW)/dlon) reports(i,jl,il) = reports(i,jl,il) + 1 end if end do end do ; Grid Hail days do i = 0,dimsizes(years)-1 do n=0,dimsizes(size_Hdays)-1 if (year_Hdays(n) .eq. years(i)) then jl = toint((lat_Hdays(n)-latS)/dlat) il = toint((lon_Hdays(n)-lonW)/dlon) Hdays(i,jl,il) = Hdays(i,jl,il) + 1 end if end do end do reports = where(reports.eq.0, reports@_FillValue,reports) Hdays = where(Hdays.eq.0, Hdays@_FillValue,Hdays) total_reports_2d = dim_sum_n_Wrap(reports,0) mean_reports_2d = dim_avg_n_Wrap(reports,0) total_Hdays_2d = dim_sum_n_Wrap(Hdays,0) mean_Hdays_2d = dim_avg_n_Wrap(Hdays,0) range_ind = ind(anom_sst.ge.range_values(0).and.anom_sst.le.range_values(1)) count_reports_2d = dim_sum_n_Wrap(reports(range_ind,:,:),0) avg_reports_2d = dim_avg_n_Wrap(reports(range_ind,:,:),0) count_Hdays_2d = dim_sum_n_Wrap(Hdays(range_ind,:,:),0) avg_Hdays_2d = dim_avg_n_Wrap(Hdays(range_ind,:,:),0) ;**Calculate anomaly anom_reports_2d = avg_reports_2d - mean_reports_2d copy_VarCoords(avg_reports_2d,anom_reports_2d) anom_Hdays_2d = avg_Hdays_2d - mean_Hdays_2d copy_VarCoords(avg_Hdays_2d,anom_Hdays_2d) ;**Calculate frequency of occurrence (how often did anomaly pattern occur during enso conditions?) temp = reports(range_ind,:,:) dims_temp = dimsizes(temp) elyears = dims_temp(0) nx = dims_temp(2) ny = dims_temp(1) freq_occ_years = new((/elyears,ny,nx/),"double") do ii = 0, elyears-1 freq_occ_years(ii,:,:) = where(anom_reports_2d.gt.0.and.temp(ii,:,:).ge.avg_reports_2d .or. anom_reports_2d.lt.0 .and. temp(ii,:,:).le.avg_reports_2d,1,freq_occ_years@_FillValue) end do freq_occ_sum = dim_sum_n_Wrap(freq_occ_years,0) freq_occ_reports = (freq_occ_sum/elyears)*100 copy_VarCoords(avg_reports_2d,freq_occ_reports) delete(temp) delete(dims_temp) delete(elyears) delete(nx) delete(ny) delete(freq_occ_years) delete(freq_occ_sum) temp = Hdays(range_ind,:,:) dims_temp = dimsizes(temp) elyears = dims_temp(0) nx = dims_temp(2) ny = dims_temp(1) freq_occ_years = new((/elyears,ny,nx/),"double") do ii = 0, elyears-1 freq_occ_years(ii,:,:) = where(anom_Hdays_2d.gt.0.and.temp(ii,:,:).ge.avg_Hdays_2d .or. anom_Hdays_2d.lt.0 .and. temp(ii,:,:).le.avg_Hdays_2d,1,freq_occ_years@_FillValue) end do freq_occ_sum = dim_sum_n_Wrap(freq_occ_years,0) freq_occ_Hdays = (freq_occ_sum/elyears)*100 copy_VarCoords(avg_Hdays_2d,freq_occ_Hdays) ;********************************* ; create plot ;******************************** plot_A=new(4,"graphic") wks = gsn_open_wks (image_type, out_fig_name ) res = True res@gsnDraw = False res@gsnFrame = False res@gsnAddCyclic = False ; plotted data are not cyclic res@mpFillOn = False ; turn off map fill res@mpCenterLonF = 265 ; Set map center res@cnFillPalette = "BlueWhiteOrangeRed" res@cnLinePalette = "BlueWhiteOrangeRed" res@cnFillOn = True ; turn on color fill res@cnLinesOn = False ; True is default res@cnLineLabelsOn = False ; True is default res@lbLabelBarOn = True ;** Add US and State borders** res@mpOutlineBoundarySets = "AllBoundaries" res@mpGeophysicalLineColor = "black" res@mpMinLatF = latS ; minimum display latitude res@mpMaxLatF = latN ; maximum display latitude res@mpMinLonF = lonW ; minimum display longitude res@mpMaxLonF = lonE ; maximum display longitude res@tmYROn = False ; Turn off tick marks on right side of figure res@tmYLOn = False ; Turn off tick marks on right side of figure res@tmXTOn = False ; Turn off tick marks on top of figure res@tmXBOn = False ; Turn off tick marks on top of figure res@cnLevelSelectionMode = "ManualLevels" res@cnMinLevelValF = -10 ; minimum contour level res@cnMaxLevelValF = 10 ; maximum contour level res@cnLevelSpacingF = 2 ; ** Frequency of occurrence contours res2 = res res2@cnFillPalette = "WhiteBlueGreenYellowRed" res2@cnLinePalette = "WhiteBlueGreenYellowRed" res2@cnFillOn = True ; turn on color fill res2@cnLinesOn = False ; True is default res2@cnLineLabelsOn = False ; True is default res2@lbLabelBarOn = True res2@cnLevelSelectionMode = "ManualLevels" res2@cnMinLevelValF = 0 ; minimum contour level res2@cnMaxLevelValF = 100 ; maximum contour level res2@cnLevelSpacingF = 10 res@txFontHeightF = 0.03 res2@txFontHeightF = 0.03 res@tiMainString = "Anomalies-Hail Reports" res2@tiMainString = "Frequency of Occurrence (%)-Hail Reports " res@gsnCenterString = "~C~ ONI Range "+range_values(0)+" to "+range_values(1) res2@gsnCenterString = "~C~ ONI Range "+range_values(0)+" to "+range_values(1) plot_A(0) = gsn_csm_contour_map_ce(wks,anom_reports_2d,res) plot_A(1) = gsn_csm_contour_map_ce(wks,freq_occ_reports,res2) res@tiMainString = "Anomalies-Hail Days" res2@tiMainString = "Frequency of Occurrence (%)-Hail Days " res@gsnCenterString = "~C~ ONI Range"+range_values(0)+" to "+range_values(1) res2@gsnCenterString = "~C~ ONI Range"+range_values(0)+" to "+range_values(1) plot_A(2) = gsn_csm_contour_map_ce(wks,anom_Hdays_2d,res) plot_A(3) = gsn_csm_contour_map_ce(wks,freq_occ_Hdays,res2) ; panel plot only resources resP = True ; modify the panel plot resP@gsnFrame = False resP@gsnMaximize = False ; resP@gsnPanelRight = 0.1 resP@lbTitleOn = True resP@lbTitleString = " " resP@gsnPanelLabelBar = False resP@txString = "Hail (> 1.0) 1980-"+yyyy+" Composites (# of years: "+elyears+")" gsn_panel(wks,plot_A,(/2,2/),resP) frame(wks) end