- 积分
- 3665
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2011-8-8
- 最后登录
- 1970-1-1
|
楼主 |
发表于 2011-12-30 17:14:57
|
显示全部楼层
;*************************************************
; 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("/home/rcp/daizhujun/wrfout_d02_2006-07-15_12:00:00.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 = 2. ; 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, \
50.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(1,:,:), 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
脚本如上
http://www.ncl.ucar.edu/Applications/Scripts/WRF_pcp_1.ncl |
|
窥视卡
雷达卡
发表于 2011-12-30 17:13:17
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
抢沙发
千斤顶
显身卡
成长值: 19710
发表于 2011-12-30 19:16:56