To whom this may concern:
When I run wrf_CrossSection4.ncl at line 121 I recieve an error. The
error says that the function wrf_contour is not defined. I am not sure
if it has something to do with the WRFUserARW.ncl that I am using. Can
anyone take a look the scripts. I using a wrf output in this process.
Thanks
Jonathan
; Example script to produce plots for a WRF real-data run,
; with the ARW coordinate dynamics option.
; Plot data on a cross section
; This script will plot data at a set angle through a specified point
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
;load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "/ptmp/jwsmith/NCLscripts/WRFUserARW.ncl"
begin
;
; The WRF ARW input file.
; This needs to have a ".nc" appended, so just do it.
dirname = "/ptmp/jwsmith/WRFChem3/run/"
filename = "wrfout_d01_2006-05-01_00:00:00.nc"
print(" reading from file "+dirname+filename)
a = addfile (dirname + filename, "r")
; We generate plots, but what kind do we prefer?
; type = "x11"
; type = "pdf"
type = "ps"
; type = "ncgm"
wks = gsn_open_wks(type,"plt_CrossSection4")
; Set some basic resources
res = True
res_at_MainTitle = "REAL-TIME WRF"
res_at_Footer = False
pltres = True
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; What times and how many time steps are in the data set?
FirstTime = True
FirstTimeMap = True
times = wrf_user_list_times(a) ; get times in the file
ntimes = dimsizes(times) ; number of times in the file
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
xlat = wrf_user_getvar(a, "XLAT",0)
xlon = wrf_user_getvar(a, "XLONG",0)
ter = wrf_user_getvar(a, "HGT",0)
do it = 0,ntimes-1,2 ; TIME LOOP
print("Working on time: " + times(it) )
res_at_TimeLabel = times(it) ; Set Valid time to use on plots
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; First get the variables we will need
tc = wrf_user_getvar(a,"tc",it) ; T in C
rh = wrf_user_getvar(a,"rh",it) ; relative humidity
z = wrf_user_getvar(a, "z",it) ; grid point height
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
do ip = 1, 3 ; we are doing 3 plots
; all with the pivot point (plane) in the center of the domain
; at angles 0, 45 and 90
; |
; |
; angle=0 is |, angle=90 is ------
; |
; |
; Build plane (pivot point) through which the cross section will go
; OR set to zero, if start and end points are specified
; IF using plane and angle, set opts in wrf_user_intrp3d to False
dimsrh = dimsizes(rh)
plane = new(2,float)
plane = (/ dimsrh(2)/2, dimsrh(1)/2 /) ; pivot point is center of domain
opts = False
if(ip .eq. 1) then
angle = 90.
X_plane = wrf_user_intrp2d(xlon,plane,angle,opts)
X_desc = "longitude"
end if
if(ip .eq. 2) then
angle = 0.
X_plane = wrf_user_intrp2d(xlat,plane,angle,opts)
X_desc = "latitude"
end if
if(ip .eq. 3) then
angle = 45.
X_plane = wrf_user_intrp2d(xlon,plane,angle,opts)
X_desc = "longitude"
end if
; X-axis lables
dimsX = dimsizes(X_plane)
xmin = X_plane(0)
xmax = X_plane(dimsX(0)-1)
xspan = dimsX(0)-1
nxlabs = floattoint( (xmax-xmin)/2 + 1)
if (FirstTimeMap) then
lat_plane = wrf_user_intrp2d(xlat,plane,angle,opts)
lon_plane = wrf_user_intrp2d(xlon,plane,angle,opts)
mpres = True
mpres_at_mpGeophysicalLineColor = "Black"
mpres_at_mpNationalLineColor = "Black"
mpres_at_mpUSStateLineColor = "Black"
mpres_at_mpGridLineColor = "Black"
mpres_at_mpLimbLineColor = "Black"
mpres_at_mpPerimLineColor = "Black"
pltres = True
pltres_at_FramePlot = False
optsM = res
optsM_at_NoHeaderFooter = True
optsM_at_cnFillOn = True
optsM_at_lbTitleOn = False
contour = wrf_contour(a,wks,ter,optsM)
plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres)
lnres = True
lnres_at_gsLineThicknessF = 3.0
lnres_at_gsLineColor = "Red"
do ii = 0,dimsX(0)-2
gsn_polyline(wks,plot,(/lon_plane(ii),lon_plane(ii+1)/),(/lat_plane(ii),lat_plane(ii+1)/),lnres)
end do
frame(wks)
delete(lon_plane)
delete(lat_plane)
pltres_at_FramePlot = True
end if
if (FirstTime) then
; THIS IS NEEDED FOR LABLES - ALWAYS DO (Z axis only needed once. X everytime plane changes)
; Y-axis labels
zmax = 6000. ; We only want to see the first 6 km
zz = wrf_user_intrp3d(z,z,"v",plane,angle,opts)
dims = dimsizes(zz)
do imax = 0,dims(0)-1
if ( .not.ismissing(zz(imax,0)) .and. zz(imax,0) .lt. zmax ) then
zmax_pos = imax
end if
end do
zspan = zmax_pos
zmin = z(0,0,0)
zmax = zz(zmax_pos,0)
zmin=zmin/1000.
zmax=zmax/1000.
nzlabs = floattoint(zmax + 1)
FirstTime = False
; END OF ALWAYS DO
end if
; Interpolate data vertically (in z)
rh_plane = wrf_user_intrp3d(rh,z,"v",plane,angle,opts)
tc_plane = wrf_user_intrp3d(tc,z,"v",plane,angle,opts)
; Options for XY Plots
opts_xy = res
opts_xy_at_tiXAxisString = X_desc
opts_xy_at_tiYAxisString = "Height (km)"
opts_xy_at_cnMissingValPerimOn = True
opts_xy_at_cnMissingValFillColor = 0
opts_xy_at_cnMissingValFillPattern = 11
opts_xy_at_tmXTOn = False
opts_xy_at_tmYROn = False
opts_xy_at_tmXBMode = "Explicit"
opts_xy_at_tmXBValues = fspan(0,xspan,nxlabs) ; Create the correct tick marks
opts_xy_at_tmXBLabels = sprintf("%.1f",fspan(xmin,xmax,nxlabs)); Create labels
opts_xy_at_tmXBLabelFontHeightF = 0.015
opts_xy_at_tmYLMode = "Explicit"
opts_xy_at_tmYLValues = fspan(0,zspan,nzlabs) ; Create the correct tick marks
opts_xy_at_tmYLLabels = sprintf("%.1f",fspan(zmin,zmax,nzlabs)); Create labels
opts_xy_at_tiXAxisFontHeightF = 0.020
opts_xy_at_tiYAxisFontHeightF = 0.020
opts_xy_at_tmXBMajorLengthF = 0.02
opts_xy_at_tmYLMajorLengthF = 0.02
opts_xy_at_tmYLLabelFontHeightF = 0.015
opts_xy_at_PlotOrientation = tc_plane_at_Orientation
; Plotting options for RH
opts_rh = opts_xy
opts_rh_at_ContourParameters = (/ 10., 90., 10. /)
opts_rh_at_pmLabelBarOrthogonalPosF = -0.1
opts_rh_at_cnFillOn = True
opts_rh_at_cnFillColors = (/"White","White","White", \
"White","Chartreuse","Green", \
"Green3","Green4", \
"ForestGreen","PaleGreen4"/)
; Plotting options for Temperature
opts_tc = opts_xy
opts_tc_at_cnInfoLabelZone = 1
opts_tc_at_cnInfoLabelSide = "Top"
opts_tc_at_cnInfoLabelPerimOn = True
opts_tc_at_cnInfoLabelOrthogonalPosF = -0.00005
opts_tc_at_ContourParameters = (/ 5. /)
; Get the contour info for the rh and temp
contour_tc = wrf_contour(a,wks,tc_plane(0:zmax_pos,:),opts_tc)
contour_rh = wrf_contour(a,wks,rh_plane(0:zmax_pos,:),opts_rh)
; MAKE PLOTS
plot = wrf_overlays(a,wks,(/contour_rh,contour_tc/),pltres)
; Delete options and fields, so we don't have carry over
delete(opts_xy)
delete(opts_tc)
delete(opts_rh)
delete(tc_plane)
delete(rh_plane)
delete(X_plane)
end do ; make next cross section
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FirstTimeMap = False
end do ; END OF TIME LOOP
end
; function wrf_user_set_xy( var:numeric, xp:numeric, yp:numeric, x1:numeric, \
; y1:numeric,angle:numeric ,opts )
; function wrf_user_intrp3d( var3d:numeric, z:numeric, plot_type:string, \
; loc_param:numeric, angle:numeric, opts:logical )
; function wrf_user_intrp2d( var2d:numeric, \
; loc_param:numeric, angle:numeric, opts:logical )
; function wrf_user_getvar( nc_file:file, variable:string, time:integer )
; function wrf_user_list_times( nc_file:file )
; function wrf_user_latlon_to_ij( nc_file:file, latitude:numeric,
; longitude:numeric )
; function wrf_contour(nc_file:file,wks[1]: graphic, data[*][*]:numeric, \
; opt_args[1]:logical)
; function wrf_vector(nc_file:file,wks[1]: graphic, data_u[*][*]:numeric, \
; data_v[*][*]:numeric, opt_args[1]:logical)
; function wrf_map(wks[1]:graphic,in_file[1]:file,opt_args[1]:logical)
; function wrf_map_overlays(in_file[1]:file,wks:graphic,base[1]:graphic, \
; plots[*]:graphic,opt_arg[1]:logical,mp_arg[1]:logical)
; function wrf_overlays(in_file[1]:file,wks:graphic, plots[*]:graphic, \
; opt_arg[1]:logical)
;
; These next three functions are obsolete as of version 5.0.1 of NCL.
; Do not use them. Use wrf_overlays and wrf_map_overlays instead.
;
; function wrf_map_zoom(wks[1]:graphic,in_file[1]:file,opt_args[1]:logical, \
; y1:integer,y2:integer,x1:integer,x2:integer)
; procedure wrf_map_overlay(wks:graphic,base[1]:graphic, \
; plots[*]:graphic, \
; opt_arg[1]:logical)
; procedure wrf_overlay(wks:graphic, plots[*]:graphic, \
; opt_arg[1]:logical)
;
; function add_white_space(str:string,maxlen:integer)
; procedure print_opts(opts_name,opts,debug)
; procedure print_header(icount:integer,debug)
;--------------------------------------------------------------------------------
;***********************************************************************;
; function : tointeger ;
; x:numeric ;
; ;
; Convert input to integer. ;
; ;
;***********************************************************************;
undef("tointeger")
function tointeger(x:numeric)
local xi
begin
if(typeof(x).eq."double")
xi = doubletointeger(x)
else
if(typeof(x).eq."float")
xi = floattointeger(x)
else
if(isatt(x,"_FillValue"))
xi = new(dimsizes(x),integer,x@_FillValue)
else
xi = new(dimsizes(x),integer)
delete(xi@_FillValue)
end if
xi = x
end if
end if
return(xi)
end
undef("wrf_user_set_xy")
function wrf_user_set_xy( var:numeric, xp:numeric, yp:numeric, x1:numeric, \
y1:numeric, angle:numeric, opts )
; mass coordinate version of ncl user routines
local dims,x,y,slope,intercept,distance,dx,dy,dxy,npts,xy
begin
; find intersection of line and domain boundaries
dims = dimsizes(var)
if (.not. opts) then ; We have a pivot point and location and
; need to calculate the start and end points of
; the cross section
if ((angle .gt. 315.) .or. (angle .lt. 45.) .or. \
((angle .gt. 135.) .and. (angle .lt. 225.)) ) then
; x = y*slope + intercept
slope = -(360.-angle)/45.
if( angle .lt. 45. ) then
slope = angle/45.
end if
if( angle .gt. 135.) then
slope = (angle-180.)/45.
end if
intercept = xp - yp*slope
; find intersections with domain boundaries
y0 = 0.
x0 = y0*slope + intercept
if( x0 .lt. 0.) then ; intersect outside of left boundary
x0 = 0.
y0 = (x0 - intercept)/slope
end if
if( x0 .gt. dims(2)-1) then ; intersect outside of right boundary
x0 = dims(2)-1
y0 = (x0 - intercept)/slope
end if
y1 = dims(1)-1. ; need to make sure this will be a float?
x1 = y1*slope + intercept
if( x1 .lt. 0.) then ; intersect outside of left boundary
x1 = 0.
y1 = (x1 - intercept)/slope
end if
if( x1 .gt. dims(2)-1) then ; intersect outside of right boundary
x1 = dims(2)-1
y1 = (x1 - intercept)/slope
end if
else
; y = x*slope + intercept
slope = (90.-angle)/45.
if( angle .gt. 225. ) then
slope = (270.-angle)/45.
end if
intercept = yp - xp*slope
; find intersections with domain boundaries
x0 = 0.
y0 = x0*slope + intercept
if( y0 .lt. 0.) then ; intersect outside of bottom boundary
y0 = 0.
x0 = (y0 - intercept)/slope
end if
if( y0 .gt. dims(1)-1) then ; intersect outside of top boundary
y0 = dims(1)-1
x0 = (y0 - intercept)/slope
end if
x1 = dims(2)-1. ; need to make sure this will be a float?
y1 = x1*slope + intercept
if( y1 .lt. 0.) then ; intersect outside of bottom boundary
y1 = 0.
x1 = (y1 - intercept)/slope
end if
if( y1 .gt. dims(1)-1) then ; intersect outside of top boundary
y1 = dims(1)-1
x1 = (y1 - intercept)/slope
end if
end if ; we have beginning and ending points
end if
if (opts) then ; We have a specified start and end point
x0 = xp
y0 = yp
if ( x1 .gt. dims(2)-1 ) then
x1 = dims(2)
end if
if ( y1 .gt. dims(1)-1 ) then
y1 = dims(1)
end if
end if
dx = x1 - x0
dy = y1 - y0
distance = (dx*dx + dy*dy)^0.5
npts = tointeger(distance)
dxy = new(1,typeof(distance))
dxy = distance/npts
xy = new((/ npts, 2 /),typeof(x1))
dx = dx/npts
dy = dy/npts
do i=0,npts-1
xy(i,0) = x0 + i*dx
xy(i,1) = y0 + i*dy
end do
; print(xy)
return(xy)
end
;--------------------------------------------------------------------------------
undef("wrf_user_intrp3d")
function wrf_user_intrp3d( var3d:numeric, z:numeric, \
plot_type:string, \
loc_param:numeric, angle:numeric, opts:logical )
; var3d - field to inerpolate (all input fields must be unstaggered
; z - interpolate to this field (either p/z)
; plot_type - interpolate horizontally "h", or vertically "v"
; loc_param - level for horizontal plots (eg. 500hPa ; 3000m - scalar),
; plane for vertical plots (2 values representing an xy point
; on the model domain through which the vertical plane will pass
; OR 4 values specifying start and end values
; angle - 0.0 for horizontal plots, and
; an angle for vertical plots - 90 represent a WE cross section
; opts Used IF opts is TRUE, else use loc_param and angle to determine crosssection
begin
if(plot_type .eq. "h" ) then ; horizontal cross section needed
var2d = wrf_interp_3d_z(var3d,z,loc_param)
if ( z(1,1,1) .gt. 500.) then
var2d_at_PlotLevelID = loc_param + " hPa"
else
var2d_at_PlotLevelID = .001*loc_param + " km"
end if
end if
if(plot_type .eq. "v" ) then ; vertical cross section needed
; set vertical cross section
if (opts) then
xy = wrf_user_set_xy( z, loc_param(0)-1, loc_param(1)-1, \ ; the -1 is for NCL dimensions
loc_param(2)-1, loc_param(3)-1, \
angle, opts )
else
xy = wrf_user_set_xy( z, loc_param(0), loc_param(1), \
0.0, 0.0, angle, opts )
end if
xp = dimsizes(xy)
; first we interp the variable, then z
var2dtmp = wrf_interp_2d_xy( var3d, xy)
var2dz = wrf_interp_2d_xy( z, xy)
; interp to constant z grid
z_max = max(var2dz)
z_min = min(var2dz)
dz = 0.01*(z_max - z_min)
nlevels = 101
z_var2d = new( (/nlevels/), typeof(var2dz))
z_var2d(0) = z_min
var2d = new( (/nlevels, xp(0)/), typeof(var2dz))
if(var2dz(0,0) .gt. var2dz(1,0) ) then ; monotonically decreasing coordinate
dz = -dz
z_var2d(0) = z_max
end if
do i=1, nlevels-1
z_var2d(i) = z_var2d(0)+i*dz
end do
do i=0,xp(0)-1
var2d(:,i) = wrf_interp_1d( var2dtmp(:,i), var2dz(:,i), z_var2d)
end do
st_x = tointeger(xy(0,0)) + 1
st_y = tointeger(xy(0,1)) + 1
ed_x = tointeger(xy(xp(0)-1,0)) + 1
ed_y = tointeger(xy(xp(0)-1,1)) + 1
if (opts) then
var2d_at_Orientation = "Cross-Sesion: (" + \
st_x + "," + st_y + ") to (" + \
ed_x + "," + ed_y + ")"
else
var2d_at_Orientation = "Cross-Sesion: (" + \
st_x + "," + st_y + ") to (" + \
ed_x + "," + ed_y + ") ; center=(" + \
loc_param(0) + "," + loc_param(1) + \
") ; angle=" + angle
end if
end if
return(var2d)
end
;--------------------------------------------------------------------------------
undef("wrf_user_intrp2d")
function wrf_user_intrp2d( var2d:numeric, \
loc_param:numeric, angle:numeric, opts:logical )
; var2d - field to inerpolate
; loc_param - plane for vertical plots (2 values representing an xy point
; on the model domain through which the vertical plane will pass
; OR 4 values specifying start and end values
; angle - 0.0 for horizontal plots, and
; an angle for vertical plots - 90 represent a WE cross section
; opts Used IF opts is TRUE, else use loc_param and angle to determine crosssection
begin
dims = dimsizes(var2d)
var2dtmp = new( (/1,dims(0),dims(1)/), typeof(var2d))
var2dtmp(0,:,:) = var2d(:,:)
; set vertical cross section
if (opts) then
xy = wrf_user_set_xy( var2dtmp, \
loc_param(0)-1, loc_param(1)-1, \ ; the -1 is for NCL dimensions
loc_param(2)-1, loc_param(3)-1, \
angle, opts )
else
xy = wrf_user_set_xy( var2dtmp, \
loc_param(0), loc_param(1), \
0.0, 0.0, angle, opts )
end if
xp = dimsizes(xy)
var1dtmp = wrf_interp_2d_xy( var2dtmp, xy)
dims = dimsizes(var1dtmp)
var = new( (/dims(1)/), typeof(var1dtmp))
var(:) = var1dtmp(0,:)
return(var)
end
undef("wrf_user_getvar")
function wrf_user_getvar( nc_file:file, variable:string, time:integer )
begin
if( (variable .eq. "uvmet") ) then
;; Calculate winds rotated to earth coord.
;; Return a 4D array where uvmet(0,:,:,:) is umet and
;; uvmet(1,:,:,:) is vmet
pii = 3.14159265
radians_per_degree = pii/180.
v = nc_file->V(time,:,:,:)
dimv = dimsizes(v)
u = nc_file->U(time,:,:,:)
dimu = dimsizes(u)
uvmet = new( (/ 2, dimu(0), dimu(1), dimv(2)/), typeof(u))
map_projection = nc_file_at_MAP_PROJ
if(map_projection .eq. 0) then ; no projection
uvmet(0,:,:,:) = 0.5*(u(:,:,:dimv(2)-2) + u(:,:,1:dimv(2)-1))
uvmet(1,:,:,:) = 0.5*(v(:,:dimv(1)-2,:) + v(:,1:dimv(1)-1,:))
uvmet_at_description = " u,v met velocity"
uvmet_at_units = "m/s"
return(uvmet)
end if
cen_lat = nc_file_at_CEN_LAT
if(isatt(nc_file,"STAND_LON")) then
cen_long = nc_file_at_STAND_LON
else
cen_long = nc_file_at_CEN_LON
end if
true_lat1 = nc_file_at_TRUELAT1
true_lat2 = nc_file_at_TRUELAT2
if(isfilevar(nc_file,"XLAT"))
latitude = nc_file->XLAT(0,:,:)
longitude = nc_file->XLONG(0,:,:)
else
latitude = nc_file->XLAT_M(0,:,:)
longitude = nc_file->XLONG_M(0,:,:)
end if
cone = 1.
if( map_projection .eq. 1) then ; Lambert Conformal mapping
if( (fabs(true_lat1 - true_lat2) .gt. 0.1) .and. \
(fabs(true_lat2 - 90. ) .gt. 0.1) ) then
cone = 10^(cos(true_lat1*radians_per_degree)) \
-10^(cos(true_lat2*radians_per_degree))
cone = cone/(10^(tan(45. -fabs(true_lat1/2.)*radians_per_degree)) - \
10^(tan(45. -fabs(true_lat2/2.)*radians_per_degree)) )
else
cone = sin(fabs(true_lat1)*radians_per_degree)
end if
end if
if(map_projection .eq. 2) then ; polar steraographic
cone = 1.
end if
if(map_projection .eq. 3) then ; Mercator
cone = 0.
end if
uvmet = wrf_uvmet( u,v, latitude, longitude, cen_long, cone )
return(uvmet)
end if
if( variable .eq. "ua" ) then
; U interpolated to mass points
U = nc_file->U(time,:,:,:)
dim = dimsizes(U)
ua = 0.5*(U(:,:,:dim(2)-2) + U(:,:,1:dim(2)-1))
ua_at_description = "u Velocity"
ua_at_units = "m/s"
return(ua)
end if
if( variable .eq. "va" ) then
; V interpolated to mass points
V = nc_file->V(time,:,:,:)
dim = dimsizes(V)
va = 0.5*(V(:,:dim(1)-2,:)+V(:,1:dim(1)-1,:))
va_at_description = "v Velocity"
va_at_units = "m/s"
return(va)
end if
if( variable .eq. "wa" ) then
; W interpolated to mass points
W = nc_file->W(time,:,:,:)
dim = dimsizes(W)
wa = 0.5*(W(0:dim(0)-2,:,:)+W(1:dim(0)-1,:,:))
wa_at_description = "Vertical Velocity w"
wa_at_units = "m/s"
return(wa)
end if
if( variable .eq. "pressure" ) then
; Full model pressure [=base pressure (PB) + pertubation pressure (P)]
P = nc_file->P(time,:,:,:)
PB = nc_file->PB(time,:,:,:)
pressure = 0.01 * ( P + PB )
pressure_at_description = "Pressure"
pressure_at_units = "hPa"
return(pressure)
end if
if( variable .eq. "z" ) then
; Height [=full geopotentail height / 9.81]
PH = nc_file->PH(time,:,:,:)
PHB = nc_file->PHB(time,:,:,:)
var = ( PH + PHB ) / 9.81
dim = dimsizes(var)
z = 0.5*(var(0:dim(0)-2,:,:)+var(1:dim(0)-1,:,:))
z_at_description = "Height"
z_at_units = "m"
return(z)
end if
if( variable .eq. "th" ) then
; Potentail Temperature is model output T + 300K
T = nc_file->T(time,:,:,:)
th = T + 300.
th_at_description = "Potential Temperature (theta) "
th_at_units = "K"
return(th)
end if
if( variable .eq. "tk" ) then
;; function wrf_tk needs theta and pressure (Pa) on input and returns temperature in K on return
T = nc_file->T(time,:,:,:)
th = T + 300.
P = nc_file->P(time,:,:,:)
PB = nc_file->PB(time,:,:,:)
p = P + PB
tk = wrf_tk( p , th )
return(tk)
end if
if( variable .eq. "tc" ) then
;; function wrf_tk needs theta and pressure (Pa) on input and returns temperature in K on return
T = nc_file->T(time,:,:,:)
th = T + 300.
P = nc_file->P(time,:,:,:)
PB = nc_file->PB(time,:,:,:)
p = P + PB
tc = wrf_tk( p , th )
tc = tc - 273.16
tc_at_units = "C" ; Overwrite return units
return(tc)
end if
if( variable .eq. "td" ) then
;; function wrf_td needs qv and pressure (Pa) on input and returns dewpoint temperature on return
QVAPOR = nc_file->QVAPOR(time,:,:,:)
QVAPOR = QVAPOR > 0.0
P = nc_file->P(time,:,:,:)
PB = nc_file->PB(time,:,:,:)
p = P + PB
td = wrf_td( p , QVAPOR )
return(td)
end if
if( variable .eq. "td2" ) then
;; function wrf_td needs qv and pressure (Pa) on input and returns dewpoint temperature on return
Q2 = nc_file->Q2(time,:,:)
Q2 = Q2 > 0.0
P = nc_file->P(time,0,:,:) ; get only the lowest level, to match Q2
PB = nc_file->PB(time,0,:,:) ; get only the lowest level, to match Q2
p = P + PB
td2 = wrf_td( p , Q2 )
td2_at_description = "2m Dewpoint Temperature" ; Overwrite return description
return(td2)
end if
if( variable .eq. "slp" ) then
;; first compute theta
;; function wrf_tk needs theta and pressure (Pa) on input
T = nc_file->T(time,:,:,:)
th = T + 300.
P = nc_file->P(time,:,:,:)
PB = nc_file->PB(time,:,:,:)
p = P + PB
tk = wrf_tk( p , th )
;; now compute sea level pressure, from qv, p (Pa), tk, z
QVAPOR = nc_file->QVAPOR(time,:,:,:)
QVAPOR = QVAPOR > 0.000
PH = nc_file->PH(time,:,:,:)
PHB = nc_file->PHB(time,:,:,:)
var = ( PH + PHB ) / 9.81
dim = dimsizes(var)
z = 0.5 * ( var(0:dim(0)-2,:,:) + var(1:dim(0)-1,:,:) )
;; get slp
slp = wrf_slp( z, tk, p, QVAPOR )
return(slp)
end if
if( variable .eq. "rh" ) then
;; first compuet tk
;; function wrf_tk needs theta and pressure (Pa) on input
T = nc_file->T(time,:,:,:)
th = T + 300.
P = nc_file->P(time,:,:,:)
PB = nc_file->PB(time,:,:,:)
p = P + PB
tk = wrf_tk( p , th )
;; now compute rh, using tk, p (Pa), QVAPOR
QVAPOR = nc_file->QVAPOR(time,:,:,:)
QVAPOR = QVAPOR > 0.000
rh = wrf_rh( QVAPOR, p, tk )
return(rh)
end if
; end of diagnostic variable list - we must want a variable already in
; the file. check variable dimensionality and pull proper time out of file
ndims = dimsizes(filevardimsizes(nc_file,variable))
if( ndims .eq. 4) then
var = nc_file->$variable$(time,:,:,:)
end if
if( ndims .eq. 3) then
var = nc_file->$variable$(time,:,:)
end if
if( ndims .eq. 2) then
var = nc_file->$variable$(time,:)
end if
if( ndims .eq. 1) then
var = nc_file->$variable$(time)
end if
return(var)
end
;--------------------------------------------------------------------------------
undef("wrf_user_list_times")
function wrf_user_list_times( nc_file:file )
local times, times_in_file, dims, i
begin
times_in_file = nc_file->Times
dims = dimsizes(times_in_file)
times = new(dims(0),string)
do i=0,dims(0)-1
times(i) = chartostring(times_in_file(i,:))
end do
times_at_description = "times in file"
print(times)
return(times)
end
;--------------------------------------------------------------------------------
undef("wrf_user_latlon_to_ij")
function wrf_user_latlon_to_ij( nc_file:file, latitude:numeric, \
longitude:numeric )
begin
if(isfilevar(nc_file,"XLAT"))
XLAT = nc_file->XLAT(0,:,:)
XLONG = nc_file->XLONG(0,:,:)
else
XLAT = nc_file->XLAT_M(0,:,:)
XLONG = nc_file->XLONG_M(0,:,:)
end if
loc = wrf_latlon_to_ij( XLAT, XLONG, latitude, longitude )
return(loc)
end
;--------------------------------------------------------------------------------
;--------------------------------------------------------------------------------
;--------------------------------------------------------------------------------
undef("delete_attrs")
procedure delete_attrs(opts:logical)
; This procedure does some cleanup by removing unneeded attributes
; so they don't get passed to other routines by accident.
begin
list_attrs = (/"MainTitle","MainTitlePos","MainTitlePosF", \
"InitTime","ValidTime","TimePos","TimePosF", \
"NoHeaderFooter","TimeLabel","LevelLabel", \
"FieldTitle","UnitLabel","NumVectors","AspectRatio", \
"SubFieldTitle","PlotOrientation","PlotLevelID", \
"ContourParameters","FontHeightF","Footer"/)
do i=0,dimsizes(list_attrs)-1
if(isatt(opts,list_attrs(i))) then
delete(opts@$list_attrs(i)$)
end if
end do
end
;--------------------------------------------------------------------------------
;--------------------------------------------------------------------------------
undef("set_cn_resources")
function set_cn_resources (data[*][*]:numeric, res:logical)
begin
opts = res
; The ContourParameters resource can either be a scalar that
; represents the contour level spacing, or it can be an array
; of three elements that represent the minimum level, the maximum
; level, and the level spacing.
;
mx = max(data)
mn = min(data)
if(mn.ne.mx.and.opts.and.isatt(opts,"ContourParameters")) then
if(dimsizes(opts_at_ContourParameters) .eq. 1) then
; Only the contour interval is specified.
nlev = tointeger((mx-mn)/opts_at_ContourParameters)+1
levels = nice_mnmxintvl(mn,mx,nlev,True)
if(levels(0) .lt. 0.) then
; Set a zero contour.
nlev = tointeger(levels(0)/opts_at_ContourParameters) - 1
levels(0) = nlev*opts_at_ContourParameters
end if
nlev = tointeger((levels(1)-levels(0))/opts_at_ContourParameters)+1
levels(1) = levels(0) + nlev*opts_at_ContourParameters
levels(2) = opts_at_ContourParameters
; Min level, max level, and level spacing are specified by user.
else
if(dimsizes(opts_at_ContourParameters) .eq. 3) then
levels = opts_at_ContourParameters
else
print("wrf_contour: Warning: illegal setting for ContourParameters attribute")
end if
end if
end if
; Contour levels
if(isvar("levels")) then
opts_at_cnLevelSelectionMode = get_res_value_keep(opts, "cnLevelSelectionMode", "ManualLevels")
opts_at_cnMinLevelValF = get_res_value_keep(opts, "cnMinLevelValF", levels(0))
opts_at_cnMaxLevelValF = get_res_value_keep(opts, "cnMaxLevelValF", levels(1))
opts_at_cnLevelSpacingF = get_res_value_keep(opts, "cnLevelSpacingF",levels(2))
delete(levels)
end if
; Set the default zero line thickness to 2, and the negative contour
; line dash pattern to 1 (0 is solid).
opts_at_gsnContourZeroLineThicknessF = get_res_value_keep(opts, "gsnContourZeroLineThicknessF",2.0)
opts_at_gsnContourNegLineDashPattern = get_res_value_keep(opts, "gsnContourNegLineDashPattern",1)
; Set resources that apply for both filled and line contour plots.
opts_at_cnFillDrawOrder = get_res_value_keep(opts,"cnFillDrawOrder", "PreDraw")
opts_at_cnLineLabelAngleF = get_res_value_keep(opts,"cnLineLabelAngleF", 0.0)
opts_at_cnLineLabelFontHeightF = get_res_value_keep(opts,"cnLineLabelFontHeightF", 0.015)
opts_at_cnInfoLabelFontHeightF = get_res_value_keep(opts,"cnInfoLabelFontHeightF", 0.015)
opts_at_cnLineLabelPerimOn = get_res_value_keep(opts,"cnLineLabelPerimOn", True)
opts_at_cnInfoLabelPerimOn = get_res_value_keep(opts,"cnInfoLabelPerimOn", False)
opts_at_cnLineLabelBackgroundColor = get_res_value_keep(opts,"cnLineLabelBackgroundColor", -1)
opts_at_cnHighLabelBackgroundColor = get_res_value_keep(opts,"cnHighLabelBackgroundColor", -1)
opts_at_cnLowLabelBackgroundColor = get_res_value_keep(opts,"cnLowLabelBackgroundColor", -1)
opts_at_cnLineColor = get_res_value_keep(opts,"cnLineColor", "Black")
opts_at_cnLineLabelFontColor = opts_at_cnLineColor
opts_at_cnLineLabelPerimColor = opts_at_cnLineColor
opts_at_cnInfoLabelFontColor = opts_at_cnLineColor
opts_at_cnHighLabelFontColor = opts_at_cnLineColor
opts_at_cnLowLabelFontColor = opts_at_cnLineColor
; Set field Title and levels if available
if(isatt(opts,"FieldTitle")) then
opts_at_cnInfoLabelString = opts_at_FieldTitle + " Contours: $CMN$ to $CMX$ by $CIU$"
else
if(isatt(data,"description")) then
opts_at_cnInfoLabelString = data_at_description + " Contours: $CMN$ to $CMX$ by $CIU$"
else
opts_at_cnInfoLabelString = " Contours: $CMN$ to $CMX$ by $CIU$"
end if
end if
return(opts)
end
;--------------------------------------------------------------------------------
;--------------------------------------------------------------------------------
undef("set_lb_resources")
function set_lb_resources (data[*][*]:numeric, res:logical)
begin
opts = res
; Somewhat convoluted way to see if a labelbar is not desired.
if(check_attr(opts,"pmTickMarkDisplayMode","Never",True).or.\
check_attr(opts,"pmTickMarkDisplayMode",-1,False).or.\
check_attr(opts,"pmTickMarkDisplayMode",0,False).or. \
check_attr(opts,"lbLabelBarOn",False,False).or.\
check_attr(opts,"lbLabelBarOn",0,False)) then
lbar_on = False
else
lbar_on = True
end if
atmp = get_res_value(opts,"lbLabelBarOn",True) ; Remove this resource
delete(atmp) ; just in case.
; Possible title for the labelbar
if(isatt(opts,"FieldTitle")) then
lb_desc = opts_at_FieldTitle
else
if(isatt(data,"description")) then
lb_desc = data_at_description
else
lb_desc = ""
end if
end if
if(isatt(opts,"UnitLabel") ) then
lb_desc = lb_desc + " (" + opts_at_UnitLabel + ")"
else
if(isatt(data,"units") .and. .not.(data_at_units.eq."")) then
lb_desc = lb_desc + " (" + data_at_units + ")"
end if
end if
if(.not.isatt(opts,"cnFillColors")) then
opts_at_gsnSpreadColors = get_res_value_keep(opts, "gsnSpreadColors", True)
end if
opts_at_cnInfoLabelOn = get_res_value_keep(opts,"cnInfoLabelOn", False)
opts_at_cnLinesOn = get_res_value_keep(opts,"cnLinesOn", False)
opts_at_cnLineLabelsOn = get_res_value_keep(opts,"cnLineLabelsOn", False)
; Labelbar resources
if(lbar_on) then
opts_at_pmLabelBarDisplayMode = get_res_value_keep(opts,"pmLabelBarDisplayMode", "Always")
opts_at_pmLabelBarSide = get_res_value_keep(opts,"pmLabelBarSide", "Bottom")
opts_at_lbAutoManage = get_res_value_keep(opts,"lbAutoManage",False)
opts_at_lbOrientation = get_res_value_keep(opts,"lbOrientation", "Horizontal")
opts_at_lbPerimOn = get_res_value_keep(opts,"lbPerimOn", False)
opts_at_lbLabelJust = get_res_value_keep(opts,"lbLabelJust", "BottomCenter")
opts_at_lbLabelAutoStride = get_res_value_keep(opts,"lbLabelAutoStride",True)
opts_at_lbBoxMinorExtentF = get_res_value_keep(opts,"lbBoxMinorExtentF", 0.13)
opts_at_lbTitleFontHeightF = get_res_value_keep(opts,"lbTitleFontHeightF", 0.015)
opts_at_lbLabelFontHeightF = get_res_value_keep(opts,"lbLabelFontHeightF", 0.015)
opts_at_pmLabelBarOrthogonalPosF = get_res_value_keep(opts,"pmLabelBarOrthogonalPosF", -0.1)
opts_at_lbTitleOn = get_res_value_keep(opts,"lbTitleOn", True)
if(lb_desc.ne."" .and. opts_at_lbTitleOn) then
opts_at_lbTitleOn = get_res_value_keep(opts,"lbTitleOn", True)
opts_at_lbTitleString = get_res_value_keep(opts,"lbTitleString", lb_desc)
opts_at_lbTitleJust = get_res_value_keep(opts,"lbTitleJust", "BottomCenter")
opts_at_lbTitleOffsetF = get_res_value_keep(opts,"lbTitleOffsetF", -0.5)
else
opts_at_lbTitleOn = False
end if
end if
return(opts)
end
;--------------------------------------------------------------------------------
;--------------------------------------------------------------------------------
undef("set_title_resources")
function set_title_resources (data[*][*]:numeric, res:logical)
begin
opts = res
; Set field Title and levels if available
if(isatt(opts,"FieldTitle")) then
SubTitles = opts_at_FieldTitle
else
if(isatt(data,"description")) then
SubTitles = data_at_description
else
SubTitles = "UnKnown"
end if
end if
if(isatt(opts,"SubFieldTitle")) then
SubTitles = SubTitles + " " + opts_at_SubFieldTitle
end if
if(isatt(opts,"UnitLabel")) then
SubTitles = SubTitles + " (" + opts_at_UnitLabel + ")"
else
if(isatt(data,"units") .and. .not.(data_at_units.eq."")) then
SubTitles = SubTitles + " (" + data_at_units + ")"
end if
end if
if (isatt(opts,"PlotLevelID")) then
SubTitles = SubTitles + " at " + opts_at_PlotLevelID
else
if (isatt(data,"PlotLevelID")) then
SubTitles = SubTitles + " at " + data_at_PlotLevelID
end if
end if
opts_at_tiMainString = SubTitles
return(opts)
end
;--------------------------------------------------------------------------------
;--------------------------------------------------------------------------------
undef("set_vc_resources")
function set_vc_resources (res:logical)
begin
opts = res
if ( isatt(opts,"vpWidthF") ) then
; num_vectors is used for vcMinDistanceF and vcRefLengthF
width = opts_at_vpWidthF
num_vectors = get_res_value(opts,"NumVectors",25.0)
opts_at_vcMinDistanceF = get_res_value_keep(opts,"vcMinDistanceF", width/num_vectors)
opts_at_vcRefLengthF = get_res_value_keep(opts,"vcRefLengthF", width/num_vectors)
else
opts_at_vcMinDistanceF = get_res_value_keep(opts,"vcMinDistanceF", 0.02)
opts_at_vcRefLengthF = get_res_value_keep(opts,"vcRefLengthF", 0.02)
end if
opts_at_vcGlyphStyle = get_res_value_keep(opts,"vcGlyphStyle", "WindBarb")
opts_at_vcWindBarbColor = get_res_value_keep(opts,"vcWindBarbColor", "Black")
opts_at_vcRefAnnoOn = get_res_value_keep(opts,"vcRefAnnoOn", False)
opts_at_vcMinFracLengthF = get_res_value_keep(opts,"vcMinFracLengthF", .2)
return(opts)
end
;--------------------------------------------------------------------------------
;--------------------------------------------------------------------------------
;procedure _SetMainTitle(nc_file:file,wks[1]:graphic,cn[1]:graphic,opts)
; This procedure checks the input data for certain attributes, and
; based on those, sets MainTitle, InitTime and ValidTime
;
; Attributes recognized by this procedure:
; MainTitle (main title - top left)
; (with Init time top right)
; TimeLabel (valid time - right under init time)
; NoHeaderFooter (switch all headers and footers off - mainly for panels)
;
; If the "NoHeaderFooter" attribute exists and is set True, then
; don't create any titles.
begin
;
opts = True
if(opts.and.isatt(opts,"NoHeaderFooter").and.opts_at_NoHeaderFooter) then
return
end if
;
; Set basic plot font
;
font_height = get_res_value_keep(opts,"FontHeightF",0.01)
;
;
; If a MainTitle attribute hasn't been set, then set to "WRF"
; Also set an Initial time
;
; MAIN Header of plot
opts = True
opts_at_MainTitle = get_res_value_keep(opts,"MainTitle", " ")
opts_at_MainTitlePos = get_res_value_keep(opts,"MainTitlePos", "Left")
opts_at_InitTime = get_res_value_keep(opts,"InitTime", True)
opts_at_ValidTime = get_res_value_keep(opts,"ValidTime", True)
opts_at_TimePos = get_res_value_keep(opts,"TimePos", "Right")
opts_at_Footer = get_res_value_keep(opts,"Footer", True)
if (opts_at_MainTitlePos .eq. "Left")
opts_at_MainTitlePos = "CenterLeft"
opts_at_MainTitlePosF = 0.0
end if
if (opts_at_MainTitlePos .eq. "Center")
opts_at_MainTitlePos = "CenterCenter"
opts_at_MainTitlePosF = 0.5
end if
if (opts_at_MainTitlePos .eq. "Right")
opts_at_MainTitlePos = "CenterRight"
opts_at_MainTitlePosF = 1.0
end if
if (opts_at_TimePos .eq. "Left")
MTOPosF = 0.30
else
MTOPosF = 0.20
end if
txt0 = create "MainPlotTilte" textItemClass wks
"txString" : opts_at_MainTitle
"txFontHeightF" : font_height*1.5
end create
anno = NhlAddAnnotation(cn,txt0)
setvalues anno
"amZone" : 3
"amSide" : "Top"
"amJust" : opts_at_MainTitlePos
"amParallelPosF" : opts_at_MainTitlePosF
"amOrthogonalPosF" : MTOPosF
"amResizeNotify" : False
end setvalues
; Time information on plot
if (opts_at_TimePos .eq. "Left")
opts_at_TimePos = "CenterLeft"
opts_at_TimePosF = 0.0
if (opts_at_MainTitlePos .eq. "CenterLeft")
MTOPosF = MTOPosF - 0.05
end if
end if
if (opts_at_TimePos .eq. "Right")
opts_at_TimePos = "CenterRight"
opts_at_TimePosF = 1.0
if (opts_at_MainTitlePos .eq. "CenterRight")
MTOPosF = MTOPosF - 0.05
end if
end if
if( isatt(nc_file,"START_DATE") .and. opts_at_InitTime ) then
InitTime = "Init: " + nc_file_at_START_DATE
txt1 = create "InitTime" textItemClass ks
"txString" : InitTime
"txFontHeightF" : font_height
end create
anno = NhlAddAnnotation(cn,txt1)
setvalues anno
"amZone" : 3
"amSide" : "Top"
"amJust" : opts_at_TimePos
"amParallelPosF" : opts_at_TimePosF
"amOrthogonalPosF" : MTOPosF
"amResizeNotify" : False
end setvalues
end if
plot_narrow = False
if((opts).and.(isatt(opts,"vpWidthF")).and.(isatt(opts,"vpHeightF"))) then
ph = opts_at_vpHeightF
pw = opts_at_vpWidthF
phw = ph/pw
if ( phw .gt. 1.8 ) then
plot_narrow = True
end if
end if
if( opts_at_ValidTime .and. isatt(opts,"TimeLabel") ) then
ValidTime = "Valid: " + opts_at_TimeLabel
MTOPosF = MTOPosF - 0.03
txt2 = create "ValidTime" textItemClass wks
"txString" : ValidTime
"txFontHeightF" : font_height
end create
anno = NhlAddAnnotation(cn,txt2)
setvalues anno
"amZone" : 3
"amSide" : "Top"
"amJust" : opts_at_TimePos
"amParallelPosF" : opts_at_TimePosF
"amOrthogonalPosF" : MTOPosF
"amResizeNotify" : False
end setvalues
end if
; Add Footer if called for
if( opts_at_Footer ) then
footer1 = nc_file_at_TITLE
dis = nc_file_at_DX / 1000.0
if(isatt(nc_file,"WEST-EAST_PATCH_END_UNSTAG"))
WE = "WEST-EAST_GRID_DIMENSION"
SN = "SOUTH-NORTH_GRID_DIMENSION"
BT = "BOTTOM-TOP_GRID_DIMENSION"
footer2 = " Phys Opt = " + nc_file_at_MP_PHYSICS + \
" ; PBL Opt = " + nc_file_at_BL_PBL_PHYSICS + \
" ; Cu Opt = " + nc_file_at_CU_PHYSICS + \
" ; WE = " + nc_file@$WE$ + \
" ; SN = " + nc_file@$SN$ + \
" ; Levels = " + nc_file@$BT$ + \
" ; Dis = " + dis + "km"
else
WE = "WEST-EAST_GRID_DIMENSION"
SN = "SOUTH-NORTH_GRID_DIMENSION"
BT = "BOTTOM-TOP_GRID_DIMENSION"
footer2 = " WE = " + nc_file@$WE$ + \
" ; SN = " + nc_file@$SN$ + \
" ; Levels = " + nc_file@$BT$ + \
" ; Dis = " + dis + "km"
end if
Footer = footer1 + "~C~" + footer2
else
Footer = " "
end if
txt3 = create "Footer" textItemClass wks
"txString" : Footer
"txFontHeightF" : font_height*.9
end create
anno = NhlAddAnnotation(cn,txt3)
setvalues anno
"amZone" : 1
; "amZone" : 7
"amJust" : "TopLeft"
"amSide" : "Bottom"
"amParallelPosF" : 0.0
"amOrthogonalPosF" : -0.55
"amResizeNotify" : False
end setvalues
; Add X-setion information if needed
if(opts.and.isatt(opts,"PlotOrientation")) then
;Xsection = "Cross-Section Orientation : " + opts_at_PlotOrientation
Xsection = opts_at_PlotOrientation
txt4 = create "Xsection" textItemClass wks
"txString" : Xsection
"txFontHeightF" : font_height*.9
end create
anno = NhlAddAnnotation(cn,txt4)
setvalues anno
"amZone" : 3
"amSide" : "Top"
"amJust" : "CenterRight"
"amParallelPosF" : 1.0
"amOrthogonalPosF" : 0.005
"amResizeNotify" : False
end setvalues
end if
end
;--------------------------------------------------------------------------------
;--------------------------------------------------------------------------------
function wrf_contour_ps(nc_file:file,wks[1]: graphic, data[*][*]:numeric, \
opt_args[1]:logical)
begin
callFrame = True
if ( isatt(opt_args,"FrameIT") ) then
if ( .not.opt_args_at_FrameIT ) then
callFrame = False
end if
delete (opt_args_at_FrameIT)
end if
lat2U = nc_file->XLAT_U(0,:,:)
lon2U = nc_file->XLONG_U(0,:,:)
opts = opt_args
opts_at_sfXArray = lon2U
opts_at_sfYArray = lat2U
opts_at_sfDataArray = data
opts_at_mpProjection = "Stereographic"
opts_at_mpEllipticalBoundary = True
opts_at_mpFillOn = False
opts_at_mpGeophysicalLineColor = get_res_value_keep(opts, "mpGeophysicalLineColor","Gray")
opts_at_mpGeophysicalLineThicknessF = get_res_value_keep(opts, "mpGeophysicalLineThicknessF",2.0)
; Set the contour resources
opts = set_cn_resources(data,opts)
opts_at_cnInfoLabelFontHeightF = 0.012
opts_at_cnLineLabelPerimOn = False
opts_at_cnInfoLabelPerimOn = True
; Find out if we are working with a contour or a shaded plot
; fill_on = False : line contour plot
; fill_on = True : filled contour plot
fill_on = get_res_value_keep(opts,"cnFillOn",False)
if(fill_on) then ; set lb resources if needed
opts_at_pmLabelBarDisplayMode = get_res_value_keep(opts,"pmLabelBarDisplayMode", "Never")
opts = set_lb_resources(data,opts)
opts_at_pmLabelBarOrthogonalPosF = 0.0
opts_at_lbTitleJust = "BottomLeft"
end if
opts_at_gsnDraw = False
opts_at_gsnFrame = False
opts_at_gsnMaximize = False
delete_attrs(opts)
cn = gsn_csm_contour_map_polar(wks,data,opts) ; Create the plot.
draw(cn)
if ( callFrame ) then
frame(wks)
end if
return (cn)
end
;--------------------------------------------------------------------------------
function wrf_vector_ps(nc_file:file,wks[1]: graphic, \
data_u[*][*]:numeric, data_v[*][*]:numeric, \
opt_args[1]:logical)
begin
callFrame = True
if ( isatt(opt_args,"FrameIT") ) then
if ( .not.opt_args_at_FrameIT ) then
callFrame = False
end if
delete (opt_args_at_FrameIT)
end if
if(isfilevar(nc_file,"XLAT"))
lat2T = nc_file->XLAT(0,:,:)
lon2T = nc_file->XLONG(0,:,:)
else
lat2T = nc_file->XLAT_M(0,:,:)
lon2T = nc_file->XLONG_M(0,:,:)
end if
opts = opt_args
opts_at_vfXArray = lon2T
opts_at_vfYArray = lat2T
opts_at_vfUDataArray = data_u
opts_at_vfVDataArray = data_v
opts_at_mpProjection = "Stereographic"
opts_at_mpEllipticalBoundary = True
opts_at_mpFillOn = False
opts_at_mpGeophysicalLineColor = get_res_value_keep(opts, "mpGeophysicalLineColor","Gray")
opts_at_mpGeophysicalLineThicknessF = get_res_value_keep(opts, "mpGeophysicalLineThicknessF",2.0)
; Set vector resources
opts = set_vc_resources(opts)
opts_at_gsnDraw = False
opts_at_gsnFrame = False
opts_at_gsnMaximize = False
delete_attrs(opts)
cn = gsn_csm_vector_map_polar(wks,data_u,data_v,opts) ; Create the plot.
draw(cn)
if ( callFrame ) then
frame(wks)
end if
return (cn)
end
;--------------------------------------------------------------------------------
;undef("wrf_contour")
function wrf_contour(nc_file:file,wks[1]: graphic, data[*][*]:numeric, \
opt_args[1]:logical)
; This function creates a contour plot and adds some titles to it.
;
; 1. Determine width to height ratio of plot.
;
; 2. First determine if this is to be a filled or line
; contour plot (fill_on)
;
; 3. If the ContourParameters attribute is set, then calculate
; the contour levels.
;
; 4. Set two resources for setting the zero contour line to
; a larger thickness, and for changing the negative contour
; lines to a dashed pattern.
;
; 5. If doing a filled contour plot, set a title for the labelbar
; based on whether a units attribute is set.
;
; 6. Make a copy of the resource list, and set some additional
; resources for filled contour plots.
;
; 7. Create the contour plot, attach the titles, and draw
; and advance the frame (if requested).
local dims
begin
opts = opt_args ; Make a copy of the resource list.
if(opts.and.isatt(opts,"mpOutlineBoundarySets")) then
delete(opts_at_mpOutlineBoundarySets)
end if
; Calculate ratio of plot width and height. Note that this doesn't
; affect the setting of gsnMaximize to True, because gsnMaximize will
; retain the aspect ratio of the plot.
if(opts.and.isatt(opts,"AspectRatio")) then
ratio = opts_at_AspectRatio
else
dims = dimsizes(data)
ratio = 1.*dims(0)/dims(1)
if(ratio .gt. 1.2) then
ratio = 1.2
end if
if(ratio .lt. 0.6667) then
ratio = 0.6667
end if
end if
if(ratio .gt. 1)
width = 0.65 * 1.0/ratio
height = 0.65
else
width = 0.85
height = 0.85 * ratio
end if
opts_at_vpWidthF = get_res_value_keep(opts,"vpWidthF", width)
opts_at_vpHeightF = get_res_value_keep(opts,"vpHeightF", height)
; Set some basic contour resources
opts = set_cn_resources(data,opts)
; Find out if we are working with a contour or a shaded plot
; fill_on = False : line contour plot
; fill_on = True : filled contour plot
fill_on = get_res_value_keep(opts,"cnFillOn",False)
if(fill_on) then ; set lb resources if needed
opts = set_lb_resources(data,opts)
atmp = get_res_value(opts,"lbLabelBarOn",True) ; Remove this resource
delete(atmp) ; just in case.
end if
; Set Title resources
opts = set_title_resources(data,opts)
; Setting gsnScale to True ensures that the tickmark lengths and labels
; will be the same size on both axes.
opts_at_gsnScale = get_res_value_keep(opts,"gsnScale", True)
; The default is not to draw the plot or advance the frame, and
; to maximize the plot in the frame.
opts_at_gsnDraw = False ; Make sure don't draw or frame or,
opts_at_gsnFrame = False ; maximize, b/c we'll do this later.
opts_at_gsnMaximize = False
opts2 = opts
delete_attrs(opts2) ; Clean up.
cn = gsn_contour(wks,data,opts2) ; Create the plot.
SetMainTitle(nc_file,wks,cn,opts) ; Set some titles
opts2_at_gsnDraw = get_res_value_keep(opts2,"gsnDraw", False)
opts2_at_gsnFrame = get_res_value_keep(opts2,"gsnFrame", False)
opts2_at_gsnMaximize = get_res_value_keep(opts2,"gsnMaximize", True)
draw_and_frame(wks,cn,opts2_at_gsnDraw,opts2_at_gsnFrame,False,opts2_at_gsnMaximize)
return(cn) ; Return
end
;--------------------------------------------------------------------------------
undef("wrf_vector")
function wrf_vector(nc_file:file,wks[1]: graphic, data_u[*][*]:numeric, \
data_v[*][*]:numeric, opt_args[1]:logical)
;
; This function creates a vector plot and adds some titles to it.
;
; 1. Determine width to height ratio of plot. Will also be use
; to calculate values for vector resources later.
;
; 2. Make a copy of the resource list, and set some additional
; resources.
;
; 3. Create the vector plot, attach the titles, and draw
; and advance the frame (if requested).
local dims
begin
opts = opt_args ; Make a copy of the resource list.
if(opts.and.isatt(opts,"mpOutlineBoundarySets")) then
delete(opts_at_mpOutlineBoundarySets)
end if
;
; The ratio is used to determine the width and height of the
; plot, and also to determine the value for the vcMinDistanceF
; resource.
;
if(opts.and.isatt(opts,"AspectRatio")) then
ratio = get_res_value(opts,"AspectRatio",0.)
else
dims = dimsizes(data_u)
ratio = 1.*dims(0)/dims(1)
if(ratio .gt. 1.2) then
ratio = 1.2
end if
if(ratio .lt. 0.6667) then
ratio = 0.6667
end if
end if
if(ratio .gt. 1)
width = 0.65/ratio
height = 0.65
else
width = 0.95
height = 0.95 * ratio
end if
opts_at_vpWidthF = get_res_value_keep(opts,"vpWidthF", width)
opts_at_vpHeightF = get_res_value_keep(opts,"vpHeightF", height)
; Set Title resources
opts = set_title_resources(data_u,opts)
; Set vector resources
opts = set_vc_resources(opts)
; Setting gsnScale to True ensures that the tickmark lengths and labels
; will be the same size on both axes.
opts_at_gsnScale = get_res_value_keep(opts,"gsnScale", True)
; The default is not to draw the plot or advance the frame, and
; to maximize the plot in the frame.
opts_at_gsnDraw = False ; Make sure don't draw or frame or,
opts_at_gsnFrame = False ; maximize, b/c we'll do this later.
opts_at_gsnMaximize = False
opts2 = opts
delete_attrs(opts2) ; Clean up.
vct = gsn_vector(wks,data_u,data_v,opts2) ; Create vector plot.
_SetMainTitle(nc_file,wks,vct,opts)
opts2_at_gsnDraw = get_res_value_keep(opts,"gsnDraw", False)
opts2_at_gsnFrame = get_res_value_keep(opts,"gsnFrame", False)
opts2_at_gsnMaximize = get_res_value_keep(opts,"gsnMaximize", True)
draw_and_frame(wks,vct,opts2_at_gsnDraw,opts2_at_gsnFrame,False, \
opts2_at_gsnMaximize)
return(vct) ; Return.
end
;--------------------------------------------------------------------------------
undef("wrf_map")
function wrf_map(wks[1]:graphic,in_file[1]:file,opt_args[1]:logical)
begin
;
; This function creates a map plot, and bases the projection on
; the MAP_PROJ attribute in the given file.
;
; 1. Make a copy of the resource list, and set some resources
; common to all map projections.
;
; 2. Determine the projection being used, and set resources based
; on that projection.
;
; 3. Create the map plot, and draw and advance the frame
; (if requested).
;
if(isatt(in_file,"MAP_PROJ"))
; if(in_file_at_MAP_PROJ .eq. 0)
; print("wrf_map: Error: attribute MAP_PROJ=0.")
; print(" Data does not have a valid map projection")
; return(new(1,graphic))
; end if
;
; There's a bug with map tickmarks that you can't set the length of
; the tickmark to 0.0 in a setvalues call. You have to do it
; during the create call. Double-check, though, b/c the bug may have
; been fixed by now.
;
opts = opt_args ; Make a copy of the resource list
opts = True
;
; Set some resources depending on what kind of map projection is
; chosen.
;
; MAP_PROJ = 1 : "LambertConformal"
; MAP_PROJ = 2 : "Stereographic"
; MAP_PROJ = 3 : "Mercator"
; MAP_PROJ = 6 : "Lat/Lon"
;
; Set some resources common to all three map projections.
;
if(any(in_file_at_MAP_PROJ.eq.(/1,2,3,6/))) then
if(isfilevar(in_file,"XLAT"))
lat = in_file->XLAT(0,:,:)
lon = in_file->XLONG(0,:,:)
else
lat = in_file->XLAT_M(0,:,:)
lon = in_file->XLONG_M(0,:,:)
end if
dims = dimsizes(lat)
;
; "LowRes" is the default that NCL uses, so you don't need to
; set it here. However, if you want a higher resolution, use
; "MediumRes". If you want higher resolution for the coastlines,
; then set it to "HighRes", but then you also need to download
; the RANGS-GSHHS database. Higher resolutions take longer to
; draw.
;
opts_at_mpDataBaseVersion = get_res_value_keep(opts, "mpDataBaseVersion","MediumRes")
;opts_at_mpOutlineBoundarySets = get_res_value_keep(opts, "mpOutlineBoundarySets", "AllBoundaries")
opts_at_mpOutlineBoundarySets = get_res_value_keep(opts, "mpOutlineBoundarySets", "GeophysicalAndUSStates")
opts_at_mpPerimLineThicknessF = get_res_value_keep(opts, "mpPerimLineThicknessF", 1.0)
opts_at_tmXBLabelFontHeightF = get_res_value_keep(opts, "tmXBLabelFontHeightF", 0.01)
opts_at_tmYLLabelFontHeightF = get_res_value_keep(opts, "tmYLLabelFontHeightF", 0.01)
;
; Select portion of the map to view.
;
;if(any(in_file_at_MAP_PROJ.eq.(/1,2,3/))) then
opts_at_mpLimitMode = get_res_value_keep(opts, "mpLimitMode","Corners")
opts_at_mpLeftCornerLatF = get_res_value_keep(opts, "mpLeftCornerLatF",lat(0,0))
opts_at_mpLeftCornerLonF = get_res_value_keep(opts, "mpLeftCornerLonF",lon(0,0))
opts_at_mpRightCornerLatF = get_res_value_keep(opts, "mpRightCornerLatF", lat(dims(0)-1,dims(1)-1))
opts_at_mpRightCornerLonF = get_res_value_keep(opts, "mpRightCornerLonF", lon(dims(0)-1,dims(1)-1))
if ( isatt(opts,"ZoomIn") .and. opts_at_ZoomIn ) then
y1 = 0
x1 = 0
y2 = dims(0)-1
x2 = dims(1)-1
if ( isatt(opts,"Ystart") ) then
y1 = opts_at_Ystart
delete(opts_at_Ystart)
end if
if ( isatt(opts,"Xstart") ) then
x1 = opts_at_Xstart
delete(opts_at_Xstart)
end if
if ( isatt(opts,"Yend") ) then
if ( opts_at_Yend .le. y2 ) then
y2 = opts_at_Yend
end if
delete(opts_at_Yend)
end if
if ( isatt(opts,"Xend") ) then
if ( opts_at_Xend .le. x2 ) then
x2 = opts_at_Xend
end if
delete(opts_at_Xend)
end if
opts_at_mpLeftCornerLatF = lat(y1,x1)
opts_at_mpLeftCornerLonF = lon(y1,x1)
opts_at_mpRightCornerLatF = lat(y2,x2)
opts_at_mpRightCornerLonF = lon(y2,x2)
delete(opts_at_ZoomIn)
end if
if ( opts_at_mpRightCornerLonF .lt. 0.0 ) then
opts_at_mpRightCornerLonF = opts_at_mpRightCornerLonF + 360.0
;print(opts_at_mpLeftCornerLonF)
;print(opts_at_mpRightCornerLonF)
end if
;end if
;
; Set some other resources for line colors and grid spacing.
;
opts_at_mpGeophysicalLineColor = get_res_value_keep(opts, "mpGeophysicalLineColor","Gray")
opts_at_mpGeophysicalLineThicknessF = get_res_value_keep(opts, "mpGeophysicalLineThicknessF",0.5)
opts_at_mpGridLineColor = get_res_value_keep(opts, "mpGridLineColor","Gray")
opts_at_mpGridLineThicknessF = get_res_value_keep(opts, "mpGridLineThicknessF",0.5)
opts_at_mpGridMaskMode = get_res_value_keep(opts, "mpGridMaskMode",3)
opts_at_mpGridSpacingF = get_res_value_keep(opts, "mpGridSpacingF",10)
opts_at_mpLimbLineColor = get_res_value_keep(opts, "mpLimbLineColor","Gray")
opts_at_mpLimbLineThicknessF = get_res_value_keep(opts, "mpLimbLineThicknessF",0.5)
opts_at_mpNationalLineColor = get_res_value_keep(opts, "mpNationalLineColor","Gray")
opts_at_mpNationalLineThicknessF = get_res_value_keep(opts, "mpNationalLineThicknessF",0.5)
opts_at_mpPerimLineColor = get_res_value_keep(opts, "mpPerimLineColor","Gray")
opts_at_mpPerimOn = get_res_value_keep(opts, "mpPerimOn",True)
opts_at_mpUSStateLineColor = get_res_value_keep(opts, "mpUSStateLineColor","Gray")
opts_at_mpUSStateLineThicknessF = get_res_value_keep(opts, "mpUSStateLineThicknessF",0.5)
opts_at_pmTickMarkDisplayMode = get_res_value_keep(opts, "pmTickMarkDisplayMode","Always")
;
; Tick mark resources
;
;opts_at_tmXBMajorLengthF = get_res_value(opts, "tmXBMajorLengthF",-0.03)
;opts_at_tmYLMajorLengthF = get_res_value(opts, "tmYLMajorLengthF",-0.03)
opts_at_tmXTOn = get_res_value(opts,"tmXTOn",False)
opts_at_tmYROn = get_res_value(opts,"tmYROn",False)
opts_at_tmYRLabelsOn = get_res_value(opts,"tmYRLabelsOn",True)
opts_at_tmXBBorderOn = get_res_value(opts,"tmXBBorderOn",True)
opts_at_tmXTBorderOn = get_res_value(opts,"tmXTBorderOn",True)
opts_at_tmYLBorderOn = get_res_value(opts,"tmYLBorderOn",True)
opts_at_tmYRBorderOn = get_res_value(opts,"tmYRBorderOn",True)
;
end if
;
; LambertConformal projection
;
if(in_file_at_MAP_PROJ .eq. 1)
projection = "LambertConformal"
opts_at_mpLambertParallel1F = get_res_value_keep(opts, "mpLambertParallel1F",in_file_at_TRUELAT1)
opts_at_mpLambertParallel2F = get_res_value_keep(opts, "mpLambertParallel2F",in_file_at_TRUELAT2)
if(isatt(in_file,"STAND_LON"))
opts_at_mpLambertMeridianF = get_res_value_keep(opts, "mpLambertMeridianF",in_file_at_STAND_LON)
else
if(isatt(in_file,"CEN_LON"))
opts_at_mpLambertMeridianF = get_res_value_keep(opts, "mpLambertMeridianF",in_file_at_CEN_LON)
else
print("ERROR: Found neither STAND_LON or CEN_LON in file")
end if
end if
end if
;
; Stereographic projection
;
if(in_file_at_MAP_PROJ .eq. 2)
projection = "Stereographic"
opts_at_mpCenterLatF = get_res_value_keep(opts, "mpCenterLatF", in_file_at_CEN_LAT)
if(isatt(in_file,"STAND_LON"))
opts_at_mpCenterLonF = get_res_value_keep(opts, "mpCenterLonF",in_file_at_STAND_LON)
else
if(isatt(in_file,"CEN_LON"))
opts_at_mpCenterLonF = get_res_value_keep(opts, "mpCenterLonF",in_file_at_CEN_LON)
else
print("ERROR: Found neither STAND_LON or CEN_LON in file")
end if
end if
end if
;
; Mercator projection
;
if(in_file_at_MAP_PROJ .eq. 3)
projection = "Mercator"
opts_at_mpCenterLatF = get_res_value_keep(opts, "mpCenterLatF", 0.0)
if(isatt(in_file,"STAND_LON"))
opts_at_mpCenterLonF = get_res_value_keep(opts, "mpCenterLonF",in_file_at_STAND_LON)
else
if(isatt(in_file,"CEN_LON"))
opts_at_mpCenterLonF = get_res_value_keep(opts, "mpCenterLonF",in_file_at_CEN_LON)
else
print("ERROR: Found neither STAND_LON or CEN_LON in file")
end if
end if
end if
;
; global WRF CylindricalEquidistant
;
if(in_file_at_MAP_PROJ .eq. 6)
projection = "CylindricalEquidistant"
opts_at_mpGridSpacingF = 45
opts_at_mpCenterLonF = get_res_value_keep(opts, "mpCenterLonF",in_file_at_CEN_LON)
if( isatt(in_file,"POLE_LAT") ) then
opts_at_mpCenterRotF = get_res_value_keep(opts, "mpCenterRotF", 90.0 - in_file_at_POLE_LAT)
end if
end if
;
; CylindricalEquidistant
;
if(in_file_at_MAP_PROJ .eq. 0)
projection = "CylindricalEquidistant"
opts_at_mpCenterLatF = get_res_value_keep(opts, "mpCenterLatF", 0.0)
if(isatt(in_file,"STAND_LON"))
opts_at_mpCenterLonF = get_res_value_keep(opts, "mpCenterLonF",in_file_at_STAND_LON)
else
if(isatt(in_file,"CEN_LON"))
opts_at_mpCenterLonF = get_res_value_keep(opts, "mpCenterLonF",in_file_at_CEN_LON)
else
print("ERROR: Found neither STAND_LON or CEN_LON in file")
end if
end if
end if
else
print("wrf_map: Error: no MAP_PROJ attribute in input file")
return(new(1,graphic))
end if
;
; The default is not to draw the plot or advance the frame, and
; to maximize the plot in the frame.
;
opts_at_gsnDraw = get_res_value_keep(opts,"gsnDraw", False)
opts_at_gsnFrame = get_res_value_keep(opts,"gsnFrame", False)
opts_at_gsnMaximize = get_res_value_keep(opts,"gsnMaximize", True)
delete_attrs(opts) ; Clean up.
mp = gsn_map(wks,projection,opts) ; Create map plot.
return(mp) ; Return.
end
;--------------------------------------------------------------------------------
undef("wrf_map_overlays")
function wrf_map_overlays(in_file[1]:file, \
wks:graphic, \
plots[*]:graphic, \
opt_arg[1]:logical, \
opt_mp[1]:logical)
; This procedure takes an array of plots and overlays them on a
; base plot - map background.
;
; It will advance the plot and cleanup, unless you set the
; PanelPlot resource to True.
;
; Attributes recognized by this procedure:
; FramePlot
; PanelPlot
; NoTitles (don't do any titles)
; CommonTitle & PlotTile is used to overwrite field titles
; CommonTitle will super-seed NoTitles
;
; If FramePlot False, then Draw the plot but do not Frame.
; In this case a user want to add to the drawing, and will
; have to advance the Frame manually in the script.
;
; If the "NoTitles" attribute exists and is set True, then
; don't create the top-left titles, and leave the main titles alone.
; This resource can be useful if you are planning to panel
; the plots.
;
; If PanelPlot is set to True, then this flags to wrf_map_overlays
; that these plots are going to be eventually paneled (likely
; by gsn_panel), and hence 1) draw and frame should not be called
; (unless gsnDraw and/or gsnFrame are explicitly set to True),
; and 2) the overlays and titles should not be removed with
; NhlRemoveOverlay and NhlRemoveAnnotation.
;
begin
; Let's make the map first
base = wrf_map(wks,in_file,opt_mp)
opts = opt_arg ; Make a copy of the resource list
no_titles = get_res_value(opts,"NoTitles",False) ; Do we want field titles?
com_title = get_res_value(opts,"CommonTitle",False) ; Do we have a common title?
if ( com_title ) then
plot_title = get_res_value(opts,"PlotTitle"," ")
no_titles = True
end if
call_draw = True
call_frame = get_res_value(opts,"FramePlot",True) ; Do we want to frame the plot?
panel_plot = get_res_value(opts,"PanelPlot",False) ; Are we paneling?
opts_at_gsnMaximize = get_res_value_keep(opts,"gsnMaximize", True)
nplots = dimsizes(plots)
; font_color = "Black"
do i=0,nplots-1
if(.not.ismissing(plots(i))) then
; class_name = NhlClassName(plots(i))
; print(class_name)
; if(class_name.eq."contourPlotClass") then
; getvalues plots(i)
; "cnFillOn" : fill_on
; "cnLineColor" : line_color
; end getvalues
; if (.not.fill_on) then
; font_color = line_color
; end if
; end if
if(.not.no_titles) then
getvalues plots(i)
"tiMainString" : SubTitle
end getvalues
if(i.eq.0) then
SubTitles = SubTitle
else
SubTitles = SubTitles + "~C~" + SubTitle
end if
end if
if(com_title .and. i .eq. nplots-1) then
getvalues plots(i)
"tiMainString" : SubTitle
end getvalues
SubTitles = plot_title
end if
setvalues plots(i)
"tfDoNDCOverlay" : True
"tiMainOn" : False
end setvalues
overlay(base,plots(i))
else
print("wrf_map_overlays: Warning: overlay plot #" + i + " is not valid.")
end if
end do
if(.not.no_titles .or. com_title) then
font_height = get_res_value_keep(opts,"FontHeightF",0.01)
txt = create "map_titles" textItemClass wks
"txString" : SubTitles
"txFontHeightF" : font_height
;"txFontColor" : font_color
end create
anno = NhlAddAnnotation(base,txt)
setvalues anno
"amZone" : 3
"amJust" : "BottomLeft"
"amSide" : "Top"
"amParallelPosF" : 0.005
"amOrthogonalPosF" : 0.03
"amResizeNotify" : False
end setvalues
base_at_map_titles = anno
end if
;
; gsnDraw and gsnFrame default to False if panel plot.
;
if(panel_plot) then
call_draw = False
call_frame= False
end if
opts_at_gsnDraw = get_res_value_keep(opts,"gsnDraw", call_draw)
opts_at_gsnFrame = get_res_value_keep(opts,"gsnFrame", call_frame)
draw_and_frame(wks,base,opts_at_gsnDraw,opts_at_gsnFrame,False, \
opts_at_gsnMaximize)
if(.not.panel_plot) then
do i=0,nplots-1
if(.not.ismissing(plots(i))) then
NhlRemoveOverlay(base,plots(i),False)
else
print("wrf_remove_map_overlays: Warning: overlay plot #" + i + " is not valid.")
print(" Nothing to remove.")
end if
end do
end if
if(.not.no_titles.and..not.panel_plot) then
if(isatt(base,"map_titles")) then
NhlRemoveAnnotation(base,base_at_map_titles)
delete(base_at_map_titles)
end if
end if
return(base)
end
;--------------------------------------------------------------------------------
undef("wrf_overlays")
function wrf_overlays(in_file[1]:file, \
wks:graphic, plots[*]:graphic, \
opt_arg[1]:logical)
; This procedure takes an array of plots and overlays them.
;
; It will advance the plot and cleanup, unless you set the
; PanelPlot resource to True.
;
; Attributes recognized by this procedure:
; FramePlot
; PanelPlot
; NoTitles (don't do any titles)
; CommonTitle & PlotTile is used to overwrite field titles
; CommonTitle will super-seed NoTitles
;
; If FramePlot False, then Draw the plot but do not Frame.
; In this case a user want to add to the drawing, and will
; have to advance the Frame manually in the script.
;
; If the "NoTitles" attribute exists and is set True, then
; don't create the top-left titles, and leave the main titles alone.
; This resource can be useful if you are planning to panel
; the plots.
;
; If PanelPlot is set to True, then this flags to wrf_overlays
; that these plots are going to be eventually paneled (likely
; by gsn_panel), and hence 1) draw and frame should not be called
; (unless gsnDraw and/or gsnFrame are explicitly set to True),
; and 2) the overlays and titles should not be removed with
; NhlRemoveOverlay and NhlRemoveAnnotation.
;
begin
opts = opt_arg ; Make a copy of the resource list.
no_titles = get_res_value(opts,"NoTitles",False) ; Do we want field titles?
com_title = get_res_value(opts,"CommonTitle",False) ; Do we have a common title?
if ( com_title ) then
plot_title = get_res_value(opts,"PlotTitle"," ")
no_titles = True
end if
call_draw = True
call_frame = get_res_value(opts,"FramePlot",True) ; Do we want to frame the plot?
panel_plot = get_res_value(opts,"PanelPlot",False) ; Are we paneling?
opts_at_gsnMaximize = get_res_value_keep(opts,"gsnMaximize", True)
nplots = dimsizes(plots)
base = plots(0)
if(.not.no_titles) then
getvalues plots(0)
"tiMainString" : SubTitle
end getvalues
SubTitles = SubTitle
setvalues plots(0)
"tfDoNDCOverlay" : True
"tiMainOn" : False
end setvalues
else
setvalues plots(0)
"tfDoNDCOverlay" : True
end setvalues
end if
if (nplots.eq.1) then
blank = create "BlankPlot" logLinPlotClass wks
;"cnConstFLabelOn" : False
end create
overlay(base,blank)
end if
do i=1,nplots-1
if(.not.ismissing(plots(i))) then
if(.not.no_titles) then
getvalues plots(i)
"tiMainString" : SubTitle
end getvalues
if(i.eq.0) then
SubTitles = SubTitle
else
SubTitles = SubTitles + "~C~" + SubTitle
end if
end if
if(com_title .and. i .eq. nplots-1) then
getvalues plots(i)
"tiMainString" : SubTitle
end getvalues
SubTitles = plot_title
end if
setvalues plots(i)
"tfDoNDCOverlay" : True
"tiMainOn" : False
end setvalues
overlay(base,plots(i))
else
print("wrf_overlays: Warning: overlay plot #" + i + " is not valid.")
end if
end do
if(.not.no_titles .or. com_title) then
font_height = get_res_value_keep(opt_arg,"FontHeightF",0.01)
txt = create "map_titles" textItemClass wks
"txString" : SubTitles
"txFontHeightF" : font_height
end create
anno = NhlAddAnnotation(base,txt)
setvalues anno
"amZone" : 3
"amJust" : "BottomLeft"
"amSide" : "Top"
"amParallelPosF" : 0.005
"amOrthogonalPosF" : 0.03
"amResizeNotify" : False
end setvalues
base_at_map_titles = anno
end if
;
; gsnDraw and gsnFrame should default to True if not a panel plot.
;
if(panel_plot) then
call_draw = False
call_frame= False
end if
opts_at_gsnDraw = get_res_value_keep(opts,"gsnDraw", call_draw)
opts_at_gsnFrame = get_res_value_keep(opts,"gsnFrame", call_frame)
opts_at_gsnMaximize = get_res_value_keep(opts,"gsnMaximize", True)
draw_and_frame(wks,base,opts_at_gsnDraw,opts_at_gsnFrame,False, \
opts_at_gsnMaximize)
if(.not.no_titles.and..not.panel_plot) then
NhlRemoveAnnotation(base,base_at_map_titles)
delete(base_at_map_titles)
end if
if(.not.panel_plot) then
if ( nplots .ge. 2 ) then
do i=1,nplots-1
if(.not.ismissing(plots(i))) then
NhlRemoveOverlay(base,plots(i),False)
else
print("wrf_remove_overlays: Warning: overlay plot #" + i + " is not valid.")
print(" Nothing to remove.")
end if
end do
end if
end if
return(base)
end
;--------------------------------------------------------------------------------
undef("wrf_map_zoom")
function wrf_map_zoom(wks[1]:graphic,in_file[1]:file,opt_args[1]:logical, \
y1:integer,y2:integer,x1:integer,x2:integer)
; As of version 5.0.1, this routine is redundant. Use the special "ZoomIn"
; resource in wrf_map to accomplish the same thing. This function is
; being kept for backwards capability. There should be no need for it
; except to run old WRF-NCL codes. Do not make any changes to it except
; possibly to fix bugs.
;
begin
print("wrf_map_zoom: Warning: This function is obsolete. Consider using")
print(" the 'ZoomIn' resource in wrf_map instead.")
if(isfilevar(in_file,"XLAT"))
lat = in_file->XLAT(0,:,:)
lon = in_file->XLONG(0,:,:)
else
lat = in_file->XLAT_M(0,:,:)
lon = in_file->XLONG_M(0,:,:)
end if
opts = opt_args ; Make a copy of the resource list
opts = True
opts_at_mpLeftCornerLatF = lat(y1,x1)
opts_at_mpLeftCornerLonF = lon(y1,x1)
opts_at_mpRightCornerLatF = lat(y2,x2)
opts_at_mpRightCornerLonF = lon(y2,x2)
mz = wrf_map(wks,in_file,opts)
return(mz)
end
;--------------------------------------------------------------------------------
undef("wrf_map_overlay")
procedure wrf_map_overlay(wks:graphic,base[1]:graphic, \
plots[*]:graphic, \
opt_arg[1]:logical)
; As of version 5.0.1, this procedure is obsolete. Use wrf_map_overlays
; instead. It is being kept for backwards capability. Do not make any
; changes to it except possibly to fix bugs.
;
; This procedure takes an array of plots and overlays them on a
; base plot - map background.
;
; It will advance the plot and cleanup, unless you set the
; PanelPlot resource to True.
;
; Attributes recognized by this procedure:
; NoTitles (don't do any titles)
; PanelPlot
;
; If the "NoTitles" attribute exists and is set True, then
; don't create the top-left titles, and leave the main titles alone.
; This resource can be useful if you are planning to panel
; the plots.
;
; If PanelPlot is set to True, then this flags to wrf_map_overlay
; that these plots are going to be eventually paneled (likely
; by gsn_panel), and hence 1) draw and frame should not be called
; (unless gsnDraw and/or gsnFrame are explicitly set to True),
; and 2) the overlays and titles should not be removed with
; NhlRemoveOverlay and NhlRemoveAnnotation.
;
begin
print("wrf_map_overlay: Warning: This procedure is obsolete. Consider" )
print(" using wrf_map_overlays instead.")
opts = opt_arg ; Make a copy of the resource list
opts = True
if(opts.and.isatt(opts,"NoTitles").and.opts_at_NoTitles) then
no_titles = True
else
no_titles = False
end if
panel_plot = get_res_value(opts,"PanelPlot",False) ; Are we paneling?
nplots = dimsizes(plots)
; font_color = "Black"
do i=0,nplots-1
if(.not.ismissing(plots(i))) then
; class_name = NhlClassName(plots(i))
; print(class_name)
; if(class_name.eq."contourPlotClass") then
; getvalues plots(i)
; "cnFillOn" : fill_on
; "cnLineColor" : line_color
; end getvalues
; if (.not.fill_on) then
; font_color = line_color
; end if
; end if
if(.not.no_titles) then
getvalues plots(i)
"tiMainString" : SubTitle
end getvalues
if(i.eq.0) then
SubTitles = SubTitle
else
SubTitles = SubTitles + "~C~" + SubTitle
end if
setvalues plots(i)
"tfDoNDCOverlay" : True
"tiMainOn" : False
end setvalues
else
setvalues plots(i)
"tfDoNDCOverlay" : True
end setvalues
end if
overlay(base,plots(i))
else
print("wrf_map_overlay: Warning: overlay plot #" + i + " is not valid.")
end if
end do
if(.not.no_titles) then
font_height = get_res_value_keep(opt_arg,"FontHeightF",0.01)
txt = create "map_titles" textItemClass wks
"txString" : SubTitles
"txFontHeightF" : font_height
;"txFontColor" : font_color
end create
anno = NhlAddAnnotation(base,txt)
setvalues anno
"amZone" : 3
"amJust" : "BottomLeft"
"amSide" : "Top"
"amParallelPosF" : 0.005
"amOrthogonalPosF" : 0.03
"amResizeNotify" : False
end setvalues
base_at_map_titles = anno
end if
;
; gsnDraw and gsnFrame should default to True if not a panel plot.
; gsnFrame will default to False if opt_arg is False.
;
if(panel_plot.or..not.opt_arg) then
call_frame = False
else
call_frame = True
end if
if(panel_plot) then
call_draw = False
else
call_draw = True
end if
opts_at_gsnDraw = get_res_value_keep(opts,"gsnDraw", call_draw)
opts_at_gsnFrame = get_res_value_keep(opts,"gsnFrame", call_frame)
opts_at_gsnMaximize = get_res_value_keep(opts,"gsnMaximize", True)
draw_and_frame(wks,base,opts_at_gsnDraw,opts_at_gsnFrame,False, \
opts_at_gsnMaximize)
if(.not.panel_plot) then
do i=0,nplots-1
if(.not.ismissing(plots(i))) then
NhlRemoveOverlay(base,plots(i),False)
else
print("wrf_remove_map_overlay: Warning: overlay plot #" + i + " is not valid.")
print(" Nothing to remove.")
end if
end do
end if
if(.not.no_titles.and..not.panel_plot) then
if(isatt(base,"map_titles")) then
NhlRemoveAnnotation(base,base_at_map_titles)
delete(base_at_map_titles)
end if
end if
end
;--------------------------------------------------------------------------------
undef("wrf_overlay")
procedure wrf_overlay(wks:graphic, plots[*]:graphic, \
opt_arg[1]:logical)
; As of version 5.0.1, this procedure is obsolete. Use wrf_overlays
; instead. It is being kept for backwards capability. Do not make any
; changes to it except possibly to fix bugs.
;
; This procedure takes an array of plots and overlays them.
;
; It will advance the plot and cleanup, unless you set the
; PanelPlot resource to True.
;
; Attributes recognized by this procedure:
; NoTitles (don't do any titles)
; PanelPlot
;
; If the "NoTitles" attribute exists and is set True, then
; don't create the top-left titles, and leave the main titles alone.
; This resource can be useful if you are planning to panel
; the plots.
;
; If PanelPlot is set to True, then this flags to wrf_overlay
; that these plots are going to be eventually paneled (likely
; by gsn_panel), and hence 1) draw and frame should not be called
; (unless gsnDraw and/or gsnFrame are explicitly set to True),
; and 2) the overlays and titles should not be removed with
; NhlRemoveOverlay and NhlRemoveAnnotation.
;
begin
print("wrf_overlay: Warning: This procedure is obsolete. Consider using")
print(" wrf_overlays instead.")
opts = opt_arg ; Make a copy of the resource list.
opts = True
if(opts.and.isatt(opts,"NoTitles").and.opts_at_NoTitles) then
no_titles = True
else
no_titles = False
end if
panel_plot = get_res_value(opts,"PanelPlot",False) ; Are we paneling?
nplots = dimsizes(plots)
base = plots(0)
if(.not.no_titles) then
getvalues plots(0)
"tiMainString" : SubTitle
end getvalues
SubTitles = SubTitle
setvalues plots(0)
"tfDoNDCOverlay" : True
"tiMainOn" : False
end setvalues
else
setvalues plots(0)
"tfDoNDCOverlay" : True
end setvalues
end if
if (nplots.eq.1) then
blank = create "BlankPlot" logLinPlotClass wks
;"cnConstFLabelOn" : False
end create
overlay(base,blank)
end if
do i=1,nplots-1
if(.not.ismissing(plots(i))) then
if(.not.no_titles) then
getvalues plots(i)
"tiMainString" : SubTitle
end getvalues
SubTitles = SubTitles + "~C~" + SubTitle
setvalues plots(i)
"tfDoNDCOverlay" : True
"tiMainOn" : False
end setvalues
else
setvalues plots(i)
"tfDoNDCOverlay" : True
end setvalues
end if
overlay(base,plots(i))
else
print("wrf_overlay: Warning: overlay plot #" + i + " is not valid.")
end if
end do
if(.not.no_titles) then
font_height = get_res_value_keep(opt_arg,"FontHeightF",0.01)
txt = create "map_titles" textItemClass wks
"txString" : SubTitles
"txFontHeightF" : font_height
end create
anno = NhlAddAnnotation(base,txt)
setvalues anno
"amZone" : 3
"amJust" : "BottomLeft"
"amSide" : "Top"
"amParallelPosF" : 0.005
"amOrthogonalPosF" : 0.03
"amResizeNotify" : False
end setvalues
base_at_map_titles = anno
end if
;
; gsnDraw and gsnFrame should default to True if not a panel plot.
; gsnFrame will default to False if opt_arg is False.
;
if(panel_plot.or..not.opt_arg) then
call_frame = False
else
call_frame = True
end if
if(panel_plot) then
call_draw = False
else
call_draw = True
end if
opts_at_gsnDraw = get_res_value_keep(opts,"gsnDraw", call_draw)
opts_at_gsnFrame = get_res_value_keep(opts,"gsnFrame", call_frame)
opts_at_gsnMaximize = get_res_value_keep(opts,"gsnMaximize", True)
draw_and_frame(wks,base,opts_at_gsnDraw,opts_at_gsnFrame,False, \
opts_at_gsnMaximize)
if(.not.no_titles.and..not.panel_plot) then
NhlRemoveAnnotation(base,base_at_map_titles)
delete(base_at_map_titles)
end if
if(.not.panel_plot) then
if ( nplots .ge. 2 ) then
do i=1,nplots-1
if(.not.ismissing(plots(i))) then
NhlRemoveOverlay(base,plots(i),False)
else
print("wrf_remove_overlay: Warning: overlay plot #" + i + " is not valid.")
print(" Nothing to remove.")
end if
end do
end if
end if
end
;--------------------------------------------------------------------------------
undef("add_white_space")
function add_white_space(str:string,maxlen:integer)
begin
cstr = stringtochar(str)
len = dimsizes(cstr)-1
ws = ""
if(len.lt.maxlen) then
do i=1,maxlen-len
ws = ws + " "
end do
end if
return(ws)
end
;--------------------------------------------------------------------------------
undef("print_opts")
procedure print_opts(opts_name,opts,debug)
begin
if(.not.debug) then
return
end if
varatts = getvaratts(opts)
;
; Sort attributes alphabetically/
;
sqsort(varatts)
;
; Get length of longest attribute name.
;
cvaratts = stringtochar(varatts)
cmaxlen = dimsizes(cvaratts(0,:))-1
print("------------------------------------------------------------")
print(opts_name + "...") ; Print name of option variable.
;
; Loop through each attribute in the list. If its value is an array,
; then print out the array with '(/' and '/)' around it.
;
; If the value is a string, then put ticks (') around each string.
;
do i=0,dimsizes(varatts)-1
x = opts@$varatts(i)$
;
; Use add_white_space to make sure all the equal signs line up.
;
tmp_str = " " + varatts(i) + \
add_white_space(varatts(i),cmaxlen) + " = "
;
; Check if attribute is an array.
;
if(dimsizes(x).gt.1) then
tmp_str = tmp_str + "(/"
do j=0,dimsizes(x)-1
if(typeof(x(j)).eq."string") then
tmp_str = tmp_str + "'" + x(j) + "'"
else
tmp_str = tmp_str + x(j)
end if
if(j.lt.dimsizes(x)-1) then
tmp_str = tmp_str + ","
else
tmp_str = tmp_str + "/)"
end if
end do
else if(typeof(x).eq."string") then
tmp_str = tmp_str + "'" + x + "'"
else
tmp_str = tmp_str + x
end if
end if
print("" + tmp_str)
delete(x)
end do
end
;--------------------------------------------------------------------------------
undef("print_header")
procedure print_header(icount:integer,debug)
begin
icount = icount + 1
if(.not.debug) then
return
end if
print("END plot #" + icount)
print("------------------------------------------------------------")
end
;--------------------------------------------------------------------------------
_______________________________________________
ncl-talk mailing list
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Fri Sep 26 2008 - 16:48:44 MDT
This archive was generated by hypermail 2.2.0 : Mon Sep 29 2008 - 13:35:43 MDT