[ncl-talk] Geopotential Height Zoomed Overlay Plot Help
Kerwyn Texeira
ktish86 at gmail.com
Tue May 3 21:12:04 MDT 2016
Hi ncl-talk,
I'm trying to plot potential height at a certain pressure level across a
terrain with winds overlayed on a zoomed in plot. I'm seeing stars. I"m
getting a long message that I don't have any idea what it means.
Error message:
warning:ContourPlotDraw: out of range coordinates encountered; standard ORV
rendering method may be unreliable;
consider setting the resource trGridType to "TriangularMesh" if
coordinates contain missing values
My Script:
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"
begin
a = addfile("./wrfout_d03_2014-01-11_22:40:00.nc","r")
it = 0
hgt = wrf_user_getvar(a, "HGT", it)
hgt at lat2d = wrf_user_getvar(a, "XLAT", it)
hgt at lon2d = wrf_user_getvar(a, "XLONG", it)
u = wrf_user_getvar(a, "ua", it)
v = wrf_user_getvar(a, "va", it)
p = wrf_user_getvar(a, "pressure", it)
qc = wrf_user_getvar(a, "QCLOUD", it)
qr = wrf_user_getvar(a, "QRAIN", it)
qs = wrf_user_getvar(a, "QSNOW", it)
qi = wrf_user_getvar(a, "QICE", it)
pre = wrf_user_getvar(a, "RAINNC", it)
z = wrf_user_getvar(a, "z", it)
u_wind = wrf_user_intrp3d(u, p, "h", 600., 0.0, False)
v_wind = wrf_user_intrp3d(v, p, "h", 600., 0.0, False)
qcl = wrf_user_intrp3d(qc, p, "h", 600., 0.0, False)
qrn = wrf_user_intrp3d(qr, p, "h", 600., 0.0, False)
qsn = wrf_user_intrp3d(qs, p, "h", 600., 0.0, False)
qice = wrf_user_intrp3d(qi, p, "h", 600., 0.0, False)
z_plane = wrf_user_intrp3d(z, p, "h", 600., 0.0, False)
qcl = qcl*1000
qrn = qrn*1000
qsn = qsn*1000
qice = qice*1000
precip = pre*0.03937
froz = qice + qsn
melt = qcl + qrn
froz at lon2d = hgt at lon2d
froz at lat2d = hgt at lat2d
melt at lon2d = hgt at lon2d
melt at lat2d = hgt at lat2d
; qcl at lon2d = hgt at lon2d
; qcl at lat2d = hgt at lat2d
; qrn at lon2d = hgt at lon2d
; qrn at lat2d = hgt at lat2d
; qsn at lon2d = hgt at lon2d
; qsn at lat2d = hgt at lat2d
; qice at lon2d = hgt at lon2d
; qice at lat2d = hgt at lat2d
spd = (u_wind*u_wind + v_wind*v_wind)^(0.5) ;m/s
u_wind = u_wind*1.94384449
v_wind = v_wind*1.94384449
u_wind at lon2d = hgt at lon2d
u_wind at lat2d = hgt at lat2d
v_wind at lon2d = hgt at lon2d
v_wind at lat2d = hgt at lat2d
; spd = spd*1.94384449
; spd at units = "Wind Speed"
; spd at units = "kts"
wks_type = "png"
wks_type at wkWidth = 2500
wks_type at wkHeight = 2500
wks = gsn_open_wks(wks_type,"geo") ; send graphics to PNG file
; gsn_define_colormap(wks,"BlAqGrYeOrRe")
gsn_define_colormap(wks,"matlab_jet")
res = True
res at gsnDraw = False ; do not draw the plot
res at gsnFrame = False ; do not advance the frame
res at cnLineLabelsOn = False ; do not use line labels
res at cnFillOn = True ; color fill
res at cnLinesOn = False ; do not draw contour lines
res at tiMainString = ""
res at pmTickMarkDisplayMode = "Always"
res at mpProjection = "CylindricalEquidistant" ;The default
res at mpDataBaseVersion = "MediumRes"
res at mpOutlineOn =True
res at lbOrientation = "Vertical"
res at tiMainOffsetYF = -0.03
res at mpFillOn = False
res at mpOutlineOn = True ; turn the map outline on
res at mpMinLatF = 37.85
res at mpMaxLatF = 38.50
res at mpMinLonF = -120.0
res at mpMaxLonF = -119.35
res at gsnLeftString = ""
res at gsnCenterString = "Geopotential Height at 600 hpa on Jan 11 at
22:40UTC"
res at gsnStringFontHeightF = 0.020
res at gsnRightString = ""
res at gsnMaximize = True
res at mpShapeMode = "FreeAspect"
res at lbTitleString = "Terrain (m)"
res at lbTitlePosition = "Right"
res at lbTitleDirection = "Across"
res at lbTitleAngleF = 90.
res at lbTitleFontHeightF = 0.015
;------wind vectors
res2 = True
res2 at gsnDraw = False
res2 at gsnFrame = False
res2 at vcWindBarbLineThicknessF= 3.0
res2 at vcRefLengthF= 0.018
res2 at vcRefMagnitudeF= 10
res2 at vcMinDistanceF = 0.05
res2 at vcGlyphStyle = "WindBarb"
res2 at gsnLeftString = ""
res2 at gsnRightString = ""
res2 at vcRefAnnoOn = False
;-------Frozen
;res3 = True
;res3 at cnLineColor = "Blue"
;res3 at gsnContourNegLineDashPattern = 2.0
;res3 at gsnContourLineThicknessesScale = 2.0
;res3 at gsnContourPosLineDashPattern = 3.0
;res3 at gsnContourZeroLineThicknessF = 0.0
;solid = gsn_csm_contour(wks,froz,res3)
;------Melted
;res4 = True
;res4 at cnLineColor = "Black"
;res4 at gsnContourLineThicknessesScale = 2.0
;res4 at cnLineLabelBackgroundColor = -1 ; transparent
;res4 at cnInfoLabelOn = False
;res4 at cnLevelSelectionMode = "ExplicitLevels"
;res4 at cnLevels = (/0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5,
0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0/)
;res4 at cnLineLabelFontHeightF = 0.01
;res4 at cnLineLabelDensityF =2
;res4 at cnLineLabelsOn = True
;liquid = gsn_csm_contour(wks,melt,res4)
;------Geopotential Height
res5 = True
res5 at cnLineColor = "Black"
res5 at gsnContourLineThicknessesScale = 3.0
geo = gsn_csm_contour(wks, z_plane, res5)
contour = gsn_csm_contour_map(wks,hgt,res)
vector = gsn_csm_vector(wks,u_wind,v_wind,res2)
overlay(contour, geo)
overlay(contour, vector)
plres = True
plres at gsLineColor = "blue"
plres at gsLineThicknessF = 3.0
lat = (/38.05, 38.05/)
lon = (/-120.10, -119.35/)
a = gsn_add_polyline(wks,contour,lon,lat,plres)
pltres = True
pltres at gsLineColor = "Black"
pltres at gsLineThicknessF = 3.0
lat1 = (/38.27, 38.27/)
lon1 = (/-120.10, -119.35/)
b = gsn_add_polyline(wks,contour,lon1,lat1,pltres)
draw(contour)
frame(wks)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
end
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20160503/6e46787b/attachment.html
More information about the ncl-talk
mailing list