;************************************************* ; WRF_pcp_1.ncl ; ; Concepts illustrated: ; - Overlaying precipitation on terrain map using wrf_xxx functions ; - Creating a contour plot with two sets of filled contours ; - Plotting WRF data ; - Creating a color map using RGB triplets ; - Merging two color maps ; - Explicitly setting contour levels ; - Using "transparent" as a contour fill color ; - Adding a title to a labelbar ; - Creating horizontal and vertical labelbars ; - Adding a vertical title to a labelbar ; 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/wrf/WRFUserARW.ncl" ;---------------------------------------------------------------------- ; This function returns an earth color palette ;---------------------------------------------------------------------- function earth_pal() local cmap begin cmap = (/ (/255, 255, 255/), \ (/ 0, 0, 0/), \ (/ 25, 25, 255/), \ (/ 20, 170, 42/), \ (/ 28, 175, 35/), \ (/ 36, 179, 29/), \ (/ 44, 184, 22/), \ (/ 51, 189, 16/), \ (/ 59, 194, 9/), \ (/ 67, 198, 2/), \ (/ 70, 200, 0/), \ (/ 71, 200, 1/), \ (/ 72, 199, 1/), \ (/ 73, 199, 2/), \ (/ 75, 198, 2/), \ (/ 76, 198, 3/), \ (/ 77, 197, 3/), \ (/ 78, 197, 4/), \ (/ 79, 197, 4/), \ (/ 80, 196, 5/), \ (/ 82, 196, 5/), \ (/ 83, 195, 6/), \ (/ 84, 195, 6/), \ (/ 85, 194, 7/), \ (/ 86, 194, 7/), \ (/ 87, 194, 8/), \ (/ 89, 193, 8/), \ (/ 90, 193, 9/), \ (/ 91, 192, 9/), \ (/ 92, 192, 10/), \ (/ 93, 191, 10/), \ (/ 94, 191, 11/), \ (/ 95, 191, 11/), \ (/ 97, 190, 12/), \ (/ 98, 190, 12/), \ (/ 99, 189, 13/), \ (/100, 189, 13/), \ (/101, 188, 14/), \ (/102, 188, 14/), \ (/104, 188, 15/), \ (/105, 187, 15/), \ (/106, 187, 16/), \ (/107, 186, 16/), \ (/108, 186, 17/), \ (/109, 185, 17/), \ (/111, 185, 18/), \ (/112, 185, 18/), \ (/113, 184, 19/), \ (/114, 184, 19/), \ (/115, 183, 20/), \ (/116, 183, 20/), \ (/118, 182, 21/), \ (/119, 182, 22/), \ (/120, 182, 22/), \ (/121, 181, 23/), \ (/122, 181, 23/), \ (/123, 180, 24/), \ (/124, 180, 24/), \ (/126, 180, 25/), \ (/127, 179, 25/), \ (/128, 179, 26/), \ (/129, 178, 26/), \ (/130, 178, 27/), \ (/131, 177, 27/), \ (/133, 177, 28/), \ (/134, 177, 28/), \ (/135, 176, 29/), \ (/136, 176, 29/), \ (/137, 175, 30/), \ (/138, 175, 30/), \ (/140, 174, 31/), \ (/141, 174, 31/), \ (/142, 174, 32/), \ (/143, 173, 32/), \ (/144, 173, 33/), \ (/145, 172, 33/), \ (/146, 172, 34/), \ (/148, 171, 34/), \ (/149, 171, 35/), \ (/150, 171, 35/), \ (/151, 170, 36/), \ (/152, 170, 36/), \ (/153, 169, 37/), \ (/155, 169, 37/), \ (/156, 168, 38/), \ (/157, 168, 38/), \ (/158, 168, 39/), \ (/159, 167, 39/), \ (/160, 167, 40/), \ (/162, 166, 40/), \ (/163, 166, 41/), \ (/164, 165, 42/), \ (/165, 165, 42/), \ (/165, 165, 43/), \ (/165, 165, 44/), \ (/166, 166, 45/), \ (/166, 166, 46/), \ (/166, 166, 47/), \ (/166, 166, 48/), \ (/166, 166, 49/), \ (/167, 167, 50/), \ (/167, 167, 51/), \ (/167, 167, 52/), \ (/167, 167, 53/), \ (/168, 168, 54/), \ (/168, 168, 55/), \ (/168, 168, 56/), \ (/168, 168, 57/), \ (/169, 169, 58/), \ (/169, 169, 59/), \ (/169, 169, 60/), \ (/169, 169, 61/), \ (/169, 169, 62/), \ (/170, 170, 63/), \ (/170, 170, 64/), \ (/170, 170, 65/), \ (/170, 170, 66/), \ (/171, 171, 67/), \ (/171, 171, 68/), \ (/171, 171, 69/), \ (/171, 171, 70/), \ (/171, 171, 71/), \ (/172, 172, 72/), \ (/172, 172, 73/), \ (/172, 172, 74/), \ (/172, 172, 75/), \ (/172, 172, 76/), \ (/173, 173, 77/), \ (/173, 173, 78/), \ (/173, 173, 79/), \ (/173, 173, 80/), \ (/174, 174, 81/), \ (/174, 174, 82/), \ (/174, 174, 83/), \ (/174, 174, 84/), \ (/175, 175, 85/), \ (/175, 175, 86/), \ (/175, 175, 87/), \ (/175, 175, 88/), \ (/175, 175, 89/), \ (/176, 176, 90/), \ (/176, 176, 91/), \ (/176, 176, 92/), \ (/176, 176, 93/), \ (/177, 177, 94/), \ (/177, 177, 95/), \ (/177, 177, 96/), \ (/177, 177, 97/), \ (/177, 177, 98/), \ (/178, 178, 99/), \ (/178, 178, 100/), \ (/178, 178, 101/), \ (/178, 178, 102/), \ (/178, 178, 103/), \ (/179, 179, 104/), \ (/179, 179, 105/), \ (/179, 179, 106/), \ (/179, 179, 107/), \ (/180, 180, 108/), \ (/180, 180, 109/), \ (/180, 180, 110/), \ (/180, 180, 111/), \ (/181, 181, 112/), \ (/181, 181, 113/), \ (/181, 181, 114/), \ (/181, 181, 115/), \ (/181, 181, 116/), \ (/182, 182, 117/), \ (/182, 182, 118/), \ (/182, 182, 119/), \ (/182, 182, 120/), \ (/183, 183, 121/), \ (/183, 183, 122/), \ (/183, 183, 123/), \ (/183, 183, 124/), \ (/183, 183, 125/), \ (/184, 184, 126/), \ (/184, 184, 127/), \ (/184, 184, 128/), \ (/184, 184, 129/), \ (/184, 184, 130/), \ (/185, 185, 131/), \ (/185, 185, 132/), \ (/185, 185, 133/), \ (/185, 185, 134/), \ (/185, 185, 135/), \ (/186, 186, 135/), \ (/186, 186, 136/), \ (/186, 186, 137/), \ (/186, 186, 138/), \ (/187, 187, 139/), \ (/187, 187, 140/), \ (/187, 187, 141/), \ (/187, 187, 142/), \ (/187, 187, 143/), \ (/188, 188, 144/), \ (/188, 188, 145/), \ (/188, 188, 146/), \ (/188, 188, 147/), \ (/188, 188, 148/), \ (/189, 189, 149/), \ (/189, 189, 150/), \ (/189, 189, 151/), \ (/189, 189, 152/), \ (/190, 190, 153/), \ (/190, 190, 154/), \ (/190, 190, 155/), \ (/190, 190, 156/), \ (/190, 190, 157/), \ (/191, 191, 158/), \ (/191, 191, 159/), \ (/191, 191, 160/), \ (/191, 191, 161/), \ (/191, 191, 162/), \ (/192, 192, 162/), \ (/192, 192, 163/), \ (/192, 192, 164/), \ (/192, 192, 165/), \ (/193, 193, 166/), \ (/193, 193, 167/), \ (/193, 193, 168/), \ (/193, 193, 169/), \ (/193, 193, 170/), \ (/194, 194, 171/), \ (/194, 194, 172/), \ (/194, 194, 173/), \ (/194, 194, 174/), \ (/194, 194, 175/), \ (/195, 195, 176/), \ (/195, 195, 177/), \ (/195, 195, 178/), \ (/195, 195, 179/), \ (/196, 196, 180/), \ (/196, 196, 181/), \ (/196, 196, 182/), \ (/196, 196, 183/), \ (/196, 196, 184/), \ (/197, 197, 185/), \ (/197, 197, 186/), \ (/197, 197, 187/), \ (/197, 197, 188/), \ (/197, 197, 189/), \ (/198, 198, 189/), \ (/198, 198, 190/), \ (/198, 198, 191/), \ (/198, 198, 192/), \ (/199, 199, 193/), \ (/199, 199, 194/), \ (/199, 199, 195/), \ (/199, 199, 196/), \ (/199, 199, 197/), \ (/200, 200, 198/), \ (/200, 200, 199/), \ (/200, 200, 200/) /)/255. return(cmap) end ;---------------------------------------------------------------------- ; This function returns a precipitation color palette ;---------------------------------------------------------------------- function precip_pal() local cmap begin cmap = (/ (/ 0, 0, 0/), \ (/242, 242, 242/), \ (/154, 192, 205/), \ (/178, 223, 238/), \ (/191, 239, 255/), \ (/ 0, 235, 235/), \ (/ 0, 163, 247/), \ (/ 0, 255, 0/), \ (/ 0, 199, 0/), \ (/ 0, 143, 0/), \ (/ 0, 63, 0/), \ (/255, 255, 0/), \ (/255, 143, 0/), \ (/255, 0, 0/), \ (/215, 0, 0/), \ (/191, 0, 0/), \ (/255, 0, 255/), \ (/155, 87, 203/), \ (/ 92, 52, 176/) /) / 255. return(cmap) end ;---------------------------------------------------------------------- ; This procedure takes the number of terrain levels and the ; number of precipitation levels, and figures out which color ; indexes to use in the earth and precipitation color palettes. ; ; It then combines a subset of each color map into one ; colormap, hopefully one that is less than 256 colors. ;---------------------------------------------------------------------- procedure set_colormap(wks,ntl,npl) local ter_cmap, pcp_cmap, nter_cols, npcp_cols, fmin, fmax, fcols, cmap begin ; Retrieve the two color maps ter_cmap = earth_pal() pcp_cmap = precip_pal() nter_cols = dimsizes(ter_cmap(:,0)) npcp_cols = dimsizes(pcp_cmap(:,0)) if( (ntl+npl+3) .gt. 256) then print("set_colormap: Error: not enough colors for the number of contour levels specified") print("No color map will be created.") return end if ; ; New color map must be (npl+1) + (ntl+1) + 2 ; (for foreground/background color) ; cmap = new((/npl+ntl+4,3/),float) cmap(0,:) = (/1.,1.,1./) cmap(1,:) = (/0.,0.,0./) ; ; Based on the terrain levels desired, calculate the indices into the ; earth color map, and put them into cmap. ; fmin = 2. fmax = nter_cols-1 fcols = fspan(fmin,fmax,ntl+1) cmap(2:ntl+2,:) = ter_cmap(floattointeger(fcols + 0.5),:) ; ; Based on the precip levels desired, calculate the indices into the ; precipitation color map, and put them into cmap. ; delete(fcols) fmin = 0. fmax = npcp_cols-1 fcols = fspan(fmin,fmax,npl+1) cmap(ntl+3:,:) = pcp_cmap(floattointeger(fcols + 0.5),:) ; Now set this color map. gsn_define_colormap(wks,cmap) end ;---------------------------------------------------------------------- ; Main code ;---------------------------------------------------------------------- begin ; ; Open file and get variables. ; a = addfile("wrfout.nc","r") ; ; Terrain ; ter = wrf_user_getvar(a,"HGT",0) ter@description = "Terrain Height" ter@units = "m" ; ; Get non-convective, convective ; Calculate total precipitation ; scale = 6. ; Artificial scale factor to see all colors rain_exp = wrf_user_getvar(a,"RAINNC",-1) / 25.4 rain_con = wrf_user_getvar(a,"RAINC", -1) / 25.4 rain_tot = (rain_exp + rain_con) * scale rain_tot@description = "Total Precipitation (inches)" ; Define contour levels here so we can determine up front ; what the merged color maps should be. precip_levels = (/ .01, 0.25, 0.50, 0.75, \ 1.00, 1.25, 1.50, 1.75, \ 2.00, 2.25, 2.50, 2.75, \ 3.00, 3.25, 3.50, 3.75, \ 4.00, 4.25, 4.50, 4.75, \ 5.00 /) ter_levels = (/ .01, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, \ 15, 16 /) * 100. nt_levs = dimsizes(ter_levels) np_levs = dimsizes(precip_levels) wks = gsn_open_wks("ps","WRF_pcp_1") ; Open graphics file. set_colormap(wks,nt_levs,np_levs) ; Define color map, which is the ; merging of two colormaps ; ; Set up resource list that will be shared between the ; two wrf_contour calls. ; res = True res@gsnDraw = False res@gsnFrame = False res@cnLevelSelectionMode = "ExplicitLevels" res@cnFillOn = True res@cnLinesOn = False res@lbLabelAutoStride = True ; ; Generate plot of terrain for background. ; ; First set up resource list specific to terrain plot. ; opts_ter = res opts_ter@cnLevels = ter_levels opts_ter@cnFillColors = ispan(2,nt_levs+2,1) ; Resources to control labelbar title, size, and location. opts_ter@lbTitleString = "Terrain" opts_ter@pmLabelBarWidthF = 0.8 opts_ter@pmLabelBarHeightF = 0.35 opts_ter@pmLabelBarOrthogonalPosF = -0.18 ; Create terrain contours (no drawing done yet). contour_ter = wrf_contour(a,wks,ter,opts_ter) ; Plotting options for Precipitation opts_r = res opts_r@cnLevels = precip_levels opts_r@cnFillColors = ispan(nt_levs+3,nt_levs+np_levs+4,1) opts_r@cnFillColors(0:1) = -1 ; make 1st two contours transparent opts_r@cnSmoothingOn = True opts_r@cnSmoothingDistanceF = .005 ; Resources to control precipitation labelbar, which will be ; vertical. opts_r@lbTitleString = "Total Precipitation" opts_r@lbTitleDirection = "Down" opts_r@lbTitleJust = "CenterRight" opts_r@lbTitlePosition = "Right" opts_r@lbTitleOffsetF = 0.07 opts_r@lbOrientation = "Vertical" opts_r@pmLabelBarSide = "Right" opts_r@pmLabelBarHeightF = 0.77 opts_r@pmLabelBarWidthF = 0.11 opts_r@pmLabelBarOrthogonalPosF = 0.03 opts_r@lbBoxMinorExtentF = 0.4 ; Total Precipitation (color fill) contour_tot = wrf_contour(a, wks, rain_tot(21,:,:), opts_r) ; ; Use the special wrf_map_overlays function to figure out ; the correct map projection and do the overlay. ; ; Set up two resource lists for wrf_map_overlays. ; pltres = True mpres = True mpres@mpGeophysicalLineColor = "Black" mpres@mpNationalLineColor = "Black" mpres@mpUSStateLineColor = "Black" mpres@mpGridLineColor = "Black" mpres@mpLimbLineColor = "Black" mpres@mpPerimLineColor = "Black" ; mpres@mpGeophysicalLineThicknessF = .5 ; mpres@mpUSStateLineThicknessF = .5 ; mpres@mpDataBaseVersion = "HighRes" ; mpres@mpDataResolution = "FinestResolution" plot = wrf_map_overlays(a,wks,(/contour_ter,contour_tot/),pltres,mpres) end