[ncl-talk] Trying to plot WRF data without using built-in WRF functions

Dennis Shea shea at ucar.edu
Mon May 11 11:27:24 MDT 2015


RE:  "This is a piece that I copy into all of my map scripts. "

I think having a personal library ['AssortedStuff.ncl'] which contains
something like the following,
is a nicer approach.

undef("wrf_setup")
function wrf_setup(a, it)
local LAT,LON, lat,lon, dimll, nlat, nlon, pltres, projres
begin
    LAT = wrf_user_getvar(a,"lat", it)
    LON = wrf_user_getvar(a,"lat", it)
    rank = dimsizes(dimsizes(LAT))

    if (rank.eq.2) then   ; it=0, it=3, ... whatever
        lat = LAT
        lon = LON
    else                        ; it = -1
        lat = LAT(0,:,:)
        lon = LON(0,:,:)
    end if

    dimll = dimsizes(lat)
    nlat  = dim(0)
    nlon = dim(1)

    projres = res
    projres at mpProjection                            = "LambertConformal"
    projres at mpRightCornerLatF                 = lat(0,0)
    projres at mpRightCornerLonF                = lon(0,0)
    projres at mpLeftCornerLatF                   = lat(nlon-1,nlat-1)
    projres at mpLeftCornerLonF                  = lon(nlon-1,nlat-1)
    projres at mpLambertParallel1F               = a at TRUELAT1
    projres at mpLambertParallel2F               = a at TRUELAT2
    projres at mpLambertMeridianF               = a at STAND_LON
    projres at mpCenterLonF                         = a at CEN_LON
    projres at mpCenterLatF                          = a at CEN_LAT
    projres at mpLimitMode                           = "Corners"
    projres at pmTickMarkDisplayMode        = "Always"
    projres at mpDataBaseVersion                  = "Ncarg4_1"
    projres at mpDataSetName                      = "Earth..4"
    projres at mpOutlineBoundarySets            = "AllBoundaries" ; all
boundaries
    projres at mpNationalLineThicknessF        =  1.0
    projres at mpGeophysicalLineThicknessF  =  1.0
    projres at mpUSStateLineThicknessF        =  2.0

   return(projres)
end

On Thu, Apr 9, 2015 at 12:31 PM, Craig Tierney - NOAA Affiliate <
craig.tierney at noaa.gov> wrote:

> Alex,
>
> Thanks for the help.  This did the trick.
>
> Craig
>
>
> On Thu, Apr 9, 2015 at 11:44 AM, Alexander Schaefer <
> alexander.schaefer at mines.sdsmt.edu> wrote:
>
>> Hi Craig,
>>
>> I do a lot of specializing of my plots with my WRF data so I have strayed
>> from using the wrf functions except getvar since it has a bunch of
>> diagnostic variables.  When I do my map projection, which is also
>> LambertConformal I set all of the following resources as such:
>>
>>     lat = wrf_user_getvar(a,"lat", it)
>>     lon = wrf_user_getvar(a,"lat", it)
>>
>>     pltres = True
>>     projres = res
>>     projres at mpProjection                            = "LambertConformal"
>>     projres at mpRightCornerLatF                 = lat(0,0)
>>     projres at mpRightCornerLonF                = lon(0,0)
>>     projres at mpLeftCornerLatF                   = lat(nlon-1,nlat-1)
>>     projres at mpLeftCornerLonF                  = lon(nlon-1,nlat-1)
>>     projres at mpLambertParallel1F               = a at TRUELAT1
>>     projres at mpLambertParallel2F               = a at TRUELAT2
>>     projres at mpLambertMeridianF               = a at STAND_LON
>>     projres at mpCenterLonF                         = a at CEN_LON
>>     projres at mpCenterLatF                          = a at CEN_LAT
>>     projres at mpLimitMode                           = "Corners"
>>     projres at pmTickMarkDisplayMode        = "Always"
>>     projres at mpDataBaseVersion                  = "Ncarg4_1"
>>     projres at mpDataSetName                      = "Earth..4"
>>     projres at mpOutlineBoundarySets            = "AllBoundaries" ; all
>> boundaries
>>     projres at mpNationalLineThicknessF        =  1.0
>>     projres at mpGeophysicalLineThicknessF  =  1.0
>>     projres at mpUSStateLineThicknessF        =  2.0
>>
>> I have altered the options so you should be able to just paste it into
>> your script.  This is a piece that I copy into all of my map scripts.
>>
>> I will then usually do other resources such as contour or vector or
>> streamlines and to keep the projection in mind instead of setting
>> "contres", contour resources equal to true, set it equal to projres and
>> then continue adding resources.
>>
>> Hope that helps.
>>
>> -Alex
>>
>> On Thu, Apr 9, 2015 at 11:07 AM, Craig Tierney - NOAA Affiliate <
>> craig.tierney at noaa.gov> wrote:
>>
>>> Hello,
>>>
>>> I am trying to plot WRF data using NCL 6.1.2 without using the built-in
>>> WRF functions.  The reason I am doing this is because what I really want to
>>> do is plot streamlines, and when I had problems doing that I tried to
>>> simplify what I was doing first.
>>>
>>> I found the examples page for plotting WRF using the gsn functions.  I
>>> started with the wrf_gsn_2.ncl example.  All I did was change the WRF input
>>> file to point to mine. The output looked like the map was
>>> CylindralEquidistant (the specified mapping) but the data are in
>>> LambertConformal, so the height field is not plotted as rectangular.  So
>>> you get the curved shape of my domain on the rectangular map where the
>>> lines of latitude are horizontal.   This isn't want I want, but the script
>>> seems to be doing what it should for my domain.
>>>
>>> What I found strange is that when I looked at the output of
>>> wrf_map_resources, that it says the Projection is LambertConformal.
>>>
>>> So I tried changing the mpProjection to LambertConformal.  When I did
>>> this I got something that looked more StereoGraphic.  The plot was centered
>>> at the north pole.
>>>
>>> When I compare the output of wrf_map_resources between this script and
>>> my script that uses wrf_map_overlays, the contents are identical.
>>>
>>> So how do I get NCL to plot my data using the GSN functions in the same
>>> projection as using wrf_map_overlays?
>>>
>>> Here is the script that I ran:
>>>
>>> ;----------------------------------------------------------------------
>>> ; wrf_gsn_2.ncl
>>> ;----------------------------------------------------------------------
>>> ; Concepts illustrated:
>>> ;   - Using gsn_csm_contour_map to plot WRF-ARW data
>>> ;----------------------------------------------------------------------
>>> ; This example is similar to wrf_gsn_1.ncl, except more plot resources
>>> ; are set to:
>>> ;   - explicitly set the contour levels
>>> ;   - change the look of the map outlines
>>> ;   - change the color map
>>> ;   - make the labelbar vertical
>>> ;----------------------------------------------------------------------;
>>> 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
>>> ;---Open WRF output file.
>>>   dir      =
>>> "/scratch2/portfolios/BMC/winds/Craig.Tierney/theia/wrf/windswrf/trunk/conusplus-wrf-6day/201308010000/nowfp-nocu/wrf/"
>>>   filename = "wrfout_d01_2013-08-04_00:00:00"
>>>   a = addfile(dir + filename + ".nc","r")
>>>
>>> ;---Read terrain height and lat/lon off file.
>>>   it        = 0     ; first time step
>>>   hgt       = wrf_user_getvar(a,"HGT",it)    ; Terrain elevation
>>>   hgt at lat2d = wrf_user_getvar(a,"XLAT",it)   ; latitude/longitude
>>>   hgt at lon2d = wrf_user_getvar(a,"XLONG",it)  ; required for plotting
>>>
>>>   wks = gsn_open_wks("png","wrf_gsn")
>>>
>>> ;---Set some basic plot options
>>>   res               = True
>>>
>>>   res at gsnMaximize   = True   ; maximize plot in frame
>>>
>>>   res at tiMainString  = filename
>>>
>>>   res at cnFillOn      = True
>>>   res at cnFillPalette = "OceanLakeLandSnow"
>>>   res at cnLinesOn     = False
>>>
>>>  res at mpProjection = "CylindricalEquidistant"
>>> ;  res at mpProjection  = "LambertConformal"
>>>
>>> ;---Zoom in on plot
>>>   res at mpMinLatF     = min(hgt at lat2d)
>>>   res at mpMaxLatF     = max(hgt at lat2d)
>>>   res at mpMinLonF     = min(hgt at lon2d)
>>>   res at mpMaxLonF     = max(hgt at lon2d)
>>>
>>> ;---Additional resources desired
>>>   res at pmTickMarkDisplayMode = "Always"   ; nicer tickmarks
>>>
>>>   res at mpDataBaseVersion     = "MediumRes"       ; better and more map
>>> outlines
>>>   res at mpDataSetName         = "Earth..4"
>>>   res at mpOutlineBoundarySets = "AllBoundaries"
>>>   res at mpOutlineOn           = True
>>>
>>>   res at lbOrientation         = "Vertical"
>>>   res at tiMainOffsetYF        = -0.03           ; Move the title down
>>>
>>> ;---Change contour levels to better match the color map being used
>>>   res at cnLevelSelectionMode = "ExplicitLevels"
>>>   res at cnLevels =
>>> (/2,100,200,400,600,800,1000,1200,1400,1600,1800,2000,2200/)
>>>
>>>
>>>   contour = gsn_csm_contour_map(wks,hgt,res)
>>>
>>> ;
>>> ; This is for debugging purposes only. It shows what map resources the
>>> ; wrf_map_overlays routine would have use, if you had called that routine
>>> ; to do the plotting. This can be useful if you are trying to reproduce
>>> ; an original WRF-ARW plot.
>>> ;
>>>   dbgres = True
>>>   dbgres = wrf_map_resources(a,dbgres)
>>>   print(dbgres)
>>> end
>>>
>>> Craig
>>>
>>>
>>> _______________________________________________
>>> ncl-talk mailing list
>>> List instructions, subscriber options, unsubscribe:
>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>
>>>
>>
>
> _______________________________________________
> ncl-talk mailing list
> 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/20150511/c3b112da/attachment.html 


More information about the ncl-talk mailing list