load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" 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/shea_util.ncl" ;************************************************ ;---------------------------------------------------------------------- ; This procedure attaches lat/lon labels to a masked lambert plot ; ; You will likely need to change lat_values and/or lon_values to ; contain the locations where you want lat/lon labels. ;---------------------------------------------------------------------- procedure add_lc_labels(wks,map,minlat,maxlat,minlon,maxlon) local lat_values, nlat, lat1_ndc, lat2_ndc, lon1_ndc, lon2_ndc,slope,txres, \ lon_values, PI, RAD_TO_DEG, dum_lft, dum_rgt, dum_bot begin PI = 3.14159 RAD_TO_DEG = 180./PI ;---Pick some "nice" values for the latitude labels. lat_values = ispan(toint(minlat),toint(maxlat),10) * 1. nlat = dimsizes(lat_values) ; ; We need to get the slope of the left and right min/max longitude lines. ; Use NDC coordinates to do this. ; lat1_ndc = new(1,float) lon1_ndc = new(1,float) lat2_ndc = new(1,float) lon2_ndc = new(1,float) datatondc(map,minlon,lat_values(0),lon1_ndc,lat1_ndc) datatondc(map,minlon,lat_values(nlat-1),lon2_ndc,lat2_ndc) slope_lft = (lat2_ndc-lat1_ndc)/(lon2_ndc-lon1_ndc) datatondc(map,maxlon,lat_values(0),lon1_ndc,lat1_ndc) datatondc(map,maxlon,lat_values(nlat-1),lon2_ndc,lat2_ndc) slope_rgt = (lat2_ndc-lat1_ndc)/(lon2_ndc-lon1_ndc) ;---Set some text resources txres = True txres@txFontHeightF = 0.03 txres@txPosXF = 0.1 ; ; Loop through lat values, and attach labels to the left and ; right edges of the masked LC plot. The labels will be ; rotated to fit the line better. ; dum_lft = new(nlat,graphic) ; Dummy array to hold attached strings. dum_rgt = new(nlat,graphic) ; Dummy array to hold attached strings. do n=0,nlat-1 ; Add extra white space to labels. lat_label_rgt = " " +lat_values(n) ;---Check if North, South, or Zero if(lat_values(n).lt.0) then lat_label_lft = lat_values(n) + "S " lat_label_rgt = lat_label_rgt + "S" end if if(lat_values(n).gt.0) then lat_label_lft = lat_values(n) + "N " lat_label_rgt = lat_label_rgt + "N" end if if(lat_values(n).eq.0) then lat_label_lft = lat_values(n) + "~S~o~N~ " end if ;---Left label txres@txAngleF = RAD_TO_DEG * atan(slope_lft) - 90 dum_lft(n) = gsn_add_text(wks,map,lat_label_lft,minlon,lat_values(n),txres) ;---Right label txres@txAngleF = RAD_TO_DEG * atan(slope_rgt) + 90 dum_rgt(n) = gsn_add_text(wks,map,lat_label_rgt,maxlon,lat_values(n),txres) end do ;---------------------------------------------------------------------- ; Now do longitude labels. These are harder because we're not ; adding them to a straight line. ; ; Loop through lon values, and attach labels to the bottom edge of the ; masked LC plot. ; delete(txres@txPosXF) txres@txPosYF = -5.0 ;---Pick some "nice" values for the longitude labels. lon_values = ispan(toint(minlon+10),toint(maxlon-10),10) * 1. nlon = dimsizes(lon_values) dum_bot = new(nlon,graphic) ; Dummy array to hold attached strings. do n=0,nlon-1 ; ; For each longitude label, we need to figure out how much to rotate ; it, so get the approximate slope at that point. ; datatondc(map,lon_values(n)-0.25,minlat,lon1_ndc,lat1_ndc) datatondc(map,lon_values(n)+0.25,minlat,lon2_ndc,lat2_ndc) slope_bot = (lat1_ndc-lat2_ndc)/(lon1_ndc-lon2_ndc) txres@txAngleF = atan(slope_bot) * RAD_TO_DEG ; ; Create longitude label. Add extra carriage returns to ; move label away from plot. ; ;---Check if East, West, or Zero lon_label_bot = " " + abs(lon_values(n)) if(lon_values(n).lt.0) then lon_label_bot = lon_label_bot + "W" end if if(lon_values(n).gt.0) then lon_label_bot = lon_label_bot + "E" end if ;---Attach to map. dum_bot(n) = gsn_add_text(wks,map,lon_label_bot,lon_values(n),minlat,txres) end do end ;---------------------------------------------------------------------- ; Main code. ;---------------------------------------------------------------------- begin ;************************************************ ; read in data1 file for binary ;************************************************ nlat = 61 mlon = 61 NMOD = 8 ; principal modes model = "mfrpspa" MODEL = "MFRPSPA" season = "DJF" season1 = "djf" pModes = "01" pModes1 = 01 nf = new((/61,61/),float) nf01 = new((/5,61,61/),float) nf4d = new((/NMOD,5,61,61/),float) mf = new((/61,61/),float) mf01 = new((/5,61,61/),float) mf4d = new((/NMOD,5,61,61/),float) kf = new((/61,61/),float) kf01 = new((/5,61,61/),float) kf4d = new((/NMOD,5,61,61/),float) var1 = "mslp" var2 = "z500" var3 = "evs" fname = "output/"+model+"_"+ season + "_" + pModes +"_conOncon" wks = gsn_open_wks("pdf",fname) ; open a pdf file path1 = model+"_"+ var1 +"_"+ season1 +".score.5x61x61" path2 = model+"_"+ var2 +"_"+ season1 +".score.5x61x61" path3 = model+"_"+ var3 +"_"+ season1 +".score.5x61x61" N=0 do nmod = 0,NMOD-1,1 do ntim = 0,4,1 nf = fbindirread(path1, N, (/nlat,mlon/), "float") nf01(ntim,:,:) = nf(:,:) nf4d(nmod,ntim,:,:)=nf01(ntim,:,:) mf = fbindirread(path2, N, (/nlat,mlon/), "float") mf01(ntim,:,:) = mf(:,:) mf4d(nmod,ntim,:,:)=mf01(ntim,:,:) kf = fbindirread(path3, N, (/nlat,mlon/), "float") kf01(ntim,:,:) = kf(:,:) kf4d(nmod,ntim,:,:)=kf01(ntim,:,:) N=N+1 end do end do ; print(nf01(1,:,:)) ; show in the screen ;************************************************ ; Create lon and lat coordinates variables ;************************************************ minlat = 10. maxlat = 40. minlon = 100. maxlon = 140. lat = new(61, float) lat = fspan( 10, 40, 61 ) lat!0 = "lat" lat@long_name = "latitude" lat@units = "degrees-north" lat&lat = lat lon = new(61,float) lon = fspan( 100, 140, 61 ) lon!0 = "lon" lon@long_name = "longitude" lon@units = "degrees-east" lon&lon = lon nlat@units = "degrees_north" mlon@units = "degrees_east" nf!0 = "lat" ; name dimensions nf!1 = "lon" nf&lat = fspan( 10, 40, 61 ) nf&lon = fspan( 100, 140, 61 ) mf!0 = "lat" ; name dimensions mf!1 = "lon" mf&lat = fspan( 10, 40, 61 ) mf&lon = fspan( 100, 140, 61 ) kf!0 = "lat" ; name dimensions kf!1 = "lon" kf&lat = fspan( 10, 40, 61 ) kf&lon = fspan( 100, 140, 61 ) ; print(nf&lat) ; print(nf&lon) ;************************************************ ; plot first contour plot ;************************************************ gsn_define_colormap(wks,"nrl_sirkes") ; choose color map ; var is mslp res = True ; plot mods desired res@mpProjection = "LambertConformal" ; choose projection res@mpFillOn = False ; turn off gray map res@mpMinLatF = minlat res@mpMaxLatF = maxlat res@mpMinLonF = minlon res@mpMaxLonF = maxlon res@mpGeophysicalLineColor = "black" ; color of cont. outlines res@mpGeophysicalLineThicknessF = 1 ; thickness of outlines res@cnInfoLabelOn = False ; turn off info label res@cnLinesOn = False ; turn off contour lines res@cnLevelSelectionMode = "ManualLevels" ; set manual contour levels res@cnMinLevelValF = -0.2 ; set min contour level res@cnMaxLevelValF = 0.2 ; set max contour level res@cnLevelSpacingF = 0.05 ; set contour spacing res@cnLineLabelsOn = False res@cnLinesOn = False ; turn off contour lines res@cnFillMode = "RasterFill" ; res@cnFillColors = (/3,4,5,7,12,12,13,15,16,17/) res@cnFillColors = (/17,16,15,13,12,12,7,5,4,3/) res@gsnDraw = False ; don't draw res@gsnFrame = False ; don't advance frame res@gsnMaximize = True ; maximize image ; res@gsnSpreadColors = True ; use full range of colors ; res@gsnSpreadColorStart = 3 ; start at color ; res@gsnSpreadColorEnd = 17 ; end at color res@sfXCStartV = minlon res@sfXCEndV = maxlon res@sfYCStartV = minlat res@sfYCEndV = maxlat res@gsnMaskLambertConformal = True ; turn on lc masking res@gsnAddCyclic = False ; res@mpGridAndLimbOn = True ; res@mpGridLatSpacingF = 10 ; res@mpGridLonSpacingF = 5 res@lbLabelBarOn = False ; turn off the label bar tiMain = (/"Day-2","Day-1","Day+0","Day+1","Day+2"/) res@tiMainOn = True res@tiMainFont = "helvetica-bold" ; set Font form res@tiMainFontHeightF = 0.02 ; title font sizei plot1 = new(6,graphic) ; create graphics array plot2 = new(6,graphic) ;************************************************ ; resource list for second data array ;************************************************ ; var is z500 res1 = True res1@sfXCStartV = minlon res1@sfXCEndV = maxlon res1@sfYCStartV = minlat res1@sfYCEndV = maxlat res1@cnLevelSelectionMode = "ManualLevels" ; set manual contour levels res1@cnMinLevelValF = -0.6 ; set min contour level res1@cnMaxLevelValF = 0.6 ; set max contour level res1@cnLevelSpacingF = 0.04 ; set contour spacing res1@cnLineColor = "Black" ; color of first contour res1@cnLineThicknessF = 2.0 ; thickness of contour res1@cnLineLabelsOn = False ; Turn on line labels. res1@cnInfoLabelOn = False ; turn off contour label res1@gsnContourNegLineDashPattern = 1 ; sets negative contours to dash pattern 1 res@tiMainString = tiMain(0) plot1(0) = gsn_csm_contour_map_overlay(wks,nf4d(pModes1-1,0,:,:),mf4d(pModes1-1,0,:,:),res,res1) res@tiMainString = tiMain(3) plot1(1) = gsn_csm_contour_map_overlay(wks,nf4d(pModes1-1,3,:,:),mf4d(pModes1-1,3,:,:),res,res1) res@tiMainString = tiMain(1) plot1(2) = gsn_csm_contour_map_overlay(wks,nf4d(pModes1-1,1,:,:),mf4d(pModes1-1,1,:,:),res,res1) res@tiMainString = tiMain(2) plot1(4) = gsn_csm_contour_map_overlay(wks,nf4d(pModes1-1,2,:,:),mf4d(pModes1-1,2,:,:),res,res1) res@tiMainString = tiMain(4) res@lbLabelBarOn = True ; turn off the label bar plot1(3) = gsn_csm_contour_map_overlay(wks,nf4d(pModes1-1,4,:,:),mf4d(pModes1-1,4,:,:),res,res1) ;---Attach latitude labels ; add_lc_labels(wks,plot1(ntim),minlat,maxlat,minlon,maxlon) ;*********************************************** ; create third contour plot ;*********************************************** ; var is evs res2 = True res2@gsnDraw = False ; don't draw res2@gsnFrame = False ; don't advance frame res2@sfXCStartV = minlon res2@sfXCEndV = maxlon res2@sfYCStartV = minlat res2@sfYCEndV = maxlat res2@cnLevelSelectionMode = "ManualLevels" ; set manual contour levels res2@cnMinLevelValF = -0.5 ; set min contour level res2@cnMaxLevelValF = 0.5 ; set max contour level res2@cnLevelSpacingF = 0.05 ; set contour spacing res2@cnLineColor = "Red" ; color of first contour res2@cnLineThicknessF = 2.0 ; thickness of contour res2@cnLineLabelsOn = False ; Turn off line labels. res2@cnInfoLabelOn = False ; turn off contour label res2@gsnContourNegLineDashPattern = 1 ; sets negative contours to dash pattern 1 plot2(0) =gsn_csm_contour(wks,kf4d(pModes1-1,0,:,:) ,res2) ; create plot plot2(1) =gsn_csm_contour(wks,kf4d(pModes1-1,3,:,:) ,res2) plot2(2) =gsn_csm_contour(wks,kf4d(pModes1-1,1,:,:) ,res2) plot2(4) =gsn_csm_contour(wks,kf4d(pModes1-1,2,:,:) ,res2) plot2(3) =gsn_csm_contour(wks,kf4d(pModes1-1,4,:,:) ,res2) do ntim = 0,4 overlay(plot1(ntim),plot2(ntim)) end do ;******************************************* ; draw panel plot with title ;******************************************* pres = True ; mod panel plot pres@gsnFrame = False ; pres@gsnPanelRight = 0.5 pres@gsnMaximize = True ; maximize plots pres@txString = MODEL+"_" + season + "_" + pModes ; add common title pres@txFont = "helvetica-bold" pres@txFontHeightF = 0.025 ; font size pres@txPosXF = 0.5 ; pres@pmLabelBarParallelPosF = 0.015 gsn_panel(wks,plot1(0:5),(/3,2/),pres) frame(wks) end