[ncl-talk] Geopotential Height Contours Overlaying Issues

Kerwyn Texeira ktish86 at gmail.com
Wed May 4 15:56:50 MDT 2016


Hi ncl-talk,

I'm trying to plot geopotential heght contours with winds overlaying
terrain zoomed in.

I have two issues;

1. The geopotential contours are not smooth.  It's plotting like it's
contours at the surface and not at 600hpa

2. Also the geopotential heights are not overlaying the terrain.  Instead,
its plotting separately.

I would like to get smoother geopotential heights at 600hpa and overlay
with winds over the terrain zoomed in.

How do I overcome these two issues? Can anyone help with these please?  I
have attached the plots.

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
  ;res at trGridType = "TrianglarMesh"
 ; res at gsnAddCyclic = False

 ;------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 cnInfoLabelOn = False
;res5 at cnLineLabelsOn = False
res5 at cnLinesOn = True
res5 at gsnLeftString = ""
res5 at gsnRightString = ""
;res5 at cnLevelSelectionMode = "ManualLevels"
;res5 at cnLineColor = "Black"
;res5 at cnLineLabelBackgroundColor = -1
;res5 at cnLineLabelPlacementMode = "constant"
;res5 at cnLineLabelFontColor = "blue"
;res5 at cnFillOpacityF = 1.0
;res5 at cnLabelDrawOrder = "PostDraw"
;res5 at cnLineDrawOrder = "PostDraw"
res5 at gsnContourLineThicknessesScale = 3.0
res5 at trGridType = "TriangularMesh"
res5 at gsnAddCyclic = False
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, vector)
overlay(contour, geo)

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/20160504/546251a6/attachment-0001.html 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: geo.000001.png
Type: image/png
Size: 479872 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20160504/546251a6/attachment-0002.png 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: geo.000002.png
Type: image/png
Size: 449284 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20160504/546251a6/attachment-0003.png 


More information about the ncl-talk mailing list