[ncl-talk] Geopotential Height Zoomed Overlay Plot Help

Mary Haley haley at ucar.edu
Wed May 4 21:06:42 MDT 2016


Kerwyn,

Please respond back to the ncl-talk email list with follow-up questions,
and not just to me personally.

It helps if you can provide a clean script, because it's hard to go through
a script that has lots of code commented out.

I took your code and tried to clean it up and use it on one of my own WRF
output files.

Are are a few issues:

[1] You are not providing lat/lon coordinate arrays for all the variables
you are plotting. You have several cases of this:

  froz at lon2d = lon2d
  froz at lat2d = lat2d

which is the right idea, but you need to do this for the variables you're
actually plotting, like "z_plane":

  z_plane at lon2d = lon2d
  z_plane at lat2d = lat2d

[2] You need to always set res at gsnAddCyclic = False. You had this in a few
places, but commented out:

 ; res at gsnAddCyclic = False


You need to uncomment these for every single plot that you are creating
that will be overlaid on a map.
[3] The "geo" plot was appearing by itself because you hadn't set:

res5 at gsnDraw = False
res5 at gsnFrame = False

And, it wasn't showing up on the "contour" plot because of the issue
mentioned in [1].

[4] You have hard-coded your min/max lat/lon resources, which I don't
normally recommend. I suggest using the min/max of your lat2d/lon2d arrays:


  res at mpMinLatF     = min(lat2d) ; 37.85

  res at mpMaxLatF     = max(lat2d) ; 38.50

  res at mpMinLonF     = min(lon2d) ; -120.0

  res at mpMaxLonF     = max(lon2d) ; -119.35

If these suggestions don't fix the problem, please post back to ncl-talk.

Thanks,

--Mary


On Wed, May 4, 2016 at 3:38 PM, Kerwyn Texeira <ktish86 at gmail.com> wrote:

> Hi ncl-talk,
>
> I not getting the errors anymore but I have two problems
>
> 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
>
>
>
> On Wed, May 4, 2016 at 7:36 AM, Mary Haley <haley at ucar.edu> wrote:
>
>> Kerwyn,
>>
>> If your lat/lon coordinate arrays have missing values, then this causes
>> issues for the contouring.  Usually, by doing what the error message
>> suggests:
>>
>> res at trGridType = "TriangularMesh"
>>
>> then what it does (I think) is throw away any data values that have
>> missing lat/lon pairs, and then it uses a triangular mesh algorithm to
>> construct the contours of the remaining data.
>>
>> However, this is WRF data, and if any of your lat2d/lon2d arrays contain
>> missing values, then there's likely an issue with your WRF output file.
>>
>> I would do this to make sure your lat/lon are okay:
>>
>>   lat2d = wrf_user_getvar(a, "XLAT", it)
>>   lon2d = wrf_user_getvar(a, "XLONG", it)
>>   printMinMax(lat2d,0)
>>   printMinMax(lon2d,0)
>>   hgt at lat2d = lat2d
>>   hgt at lon2d = lon2d
>>   . . .
>>
>> If the lat/lon arrays look okay, then the next thing I would look for is
>> to do a printVarSummary on every single variable you are trying to plot, to
>> make sure it has the correct lat2d/lon2d attached. Also, for WRF data, you
>> want to set res at gsnAddCyclic = False for any data being overlaid on a
>> map, so it doesn't try to add a longitude cyclic point.
>>
>> --Mary
>>
>>
>> On Tue, May 3, 2016 at 9:12 PM, Kerwyn Texeira <ktish86 at gmail.com> wrote:
>>
>>> 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
>>>
>>>
>>>
>>> _______________________________________________
>>> ncl-talk mailing list
>>> ncl-talk at ucar.edu
>>> List instructions, subscriber options, unsubscribe:
>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>
>>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20160504/58ce5a87/attachment.html 


More information about the ncl-talk mailing list