[ncl-talk] Overlaying Contour on WRF Domain

Mary Haley haley at ucar.edu
Thu Sep 8 10:19:27 MDT 2016


Sean,

First, I recommend using "wrf_map_resources" to set the necessary WRF map
resources for you.  This will allow you to get rid of all this code:

mpres at mpProjection                 = "LambertConformal"
mpres at mpLimitMode                  = "Corners"
mpres at mpLeftCornerLatF             = lat(0,0)
mpres at mpLeftCornerLonF             = lon(0,0)
mpres at mpRightCornerLatF            = lat(nlat-1,nlon-1)
mpres at mpRightCornerLonF            = lon(nlat-1,nlon-1)
mpres at mpLambertParallel1F          = b at TRUELAT1
mpres at mpLambertParallel2F          = b at TRUELAT2
mpres at mpLambertMeridianF           = b at STAND_LON
mpres at mpDataBaseVersion            = "MediumRes"

and replace it with one line:

mpres = wrf_map_resources(b,mpres)

The other resources you are setting:

mpres at mpOutlineBoundarySets        = "GeophysicalAndUSStates"
mpres at mpFillOn                     = False
mpres at mpOutlineOn                  = True
mpres at mpCountyLineThicknessF       = 2.5
mpres at mpGeophysicalLineThicknessF  = 2.5
mpres at mpNationalLineThicknessF     = 2.5
mpres at mpProvincialLineThicknessF   = 2.5
mpres at mpUSStateLineThicknessF      = 2.5
mpres at mpCountyLineColor            = "Black"
mpres at mpGeophysicalLineColor       = "Black"
mpres at mpNationalLineColor          = "Black"
mpres at mpUSStateLineColor           = "Black"

are also set by "wrf_map_resources", but since you are using different
values than what wrf_map_resources uses, I recommend leaving these alone.
Just make sure you set all these resources *after* you call
wrf_map_resources.  You can always do a "print(mpres)" after you call
wrf_map_resources to see exactly what values it is setting for you.

As for the "ash" array you are trying to plot: please look at your latN and
lonN arrays to make sure they look correct, and that they in "degrees"
 units. That is, the lat values should range from -90 to 90, and the
longitude values from -180 to 180 or 0 to 360:

  printMinMax(latN,0)
  printMinMax(lonN,0)

You have the sfXArray and sfYArray arrays set correctly, but you do NOT
want to set this:

cnres at sfDataArray           = ash

You also do not need this, because the gsn_csm_contour function will set it
for you:

cnres at sfMissingValueV       = ash at _FillValue

I'm not sure if this is a typo:

ash_contour = gsn_csm_contour(wks,ashNN,cnres)

but shouldn't "asnNN" be "ash"?

Finally, you want to be sure that "ash" doesn't already have some kind of
coordinate arrays attached to it.

If you continue to have problems, it would really help if we could get the
latest version of your script and your data, or, if you can't provide the
data, then insert a bunch of "printVarSummary" and "printMinMax" statements
so we can see what your data looks like:

printVarSummary(ash)
printVarSummary(latN)
printVarSummary(lonN)
printMinMax(ash,0)
printMinMax(latN,0)
printMinMax(lonN,0)

On Thu, Sep 8, 2016 at 9:57 AM, Genito Amos Maure <genito.maure at uem.mz>
wrote:

> Hi Sean.
> Is that "i" before the zero in the expression
>  lat = b->XLAT(i0,:,:) ;Not a nested domain
> a typo?
>
>
>
> On Thu, 8 Sep 2016 07:48:52 -0800, Sean Egan wrote
> > Hi,
> >
> > I have an HDF file that contains volcanic ash column densities in a
> 900x3712
> array. Most of the elements are fill values. In addition, there are two
> arrays
> of latitudes and longitudes, both 900x3712.
> >
> > I am trying to use gsn_contour or gsn_csm_contour to plot the ash
> concentrations over a WRF model domain. Currently, I'm unable to get NCL to
> reference the ash array with the latitude and longitude arrays. I have
> posted
> my script below and can upload the HDF file if needed.
> >
> > Any tips as to why I cannot get NCL to reference the lat and lon arrays
> correctly would be greatly appreciated. I believe it may have something to
> do
> with the sfXArray and sfYArray resources....
> >
> > Thanks,
> >
> > Sean
> >
> > Script Start
> >
> > 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"
> > load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
> >
> > begin
> >
> > a = addfile("file1.hdf","r")
> > b = addfile("wrfout_d01_2010-04-16_00:00:00","r")
> > wks=gsn_open_wks("x11","testing")
> >
> > latS = a->pixel_latitude
> > latN = latS * 0.002746665852839747
> > lonS = a->pixel_longitude
> > lonN = lonS * 0.005493331705679495
> > ash  = a->volcld_ret_14_15_16_ash_mass_loading
> >
> > lat = b->XLAT(i0,:,:) ;Not a nested domain
> > lon = b->XLONG(0,:,:)
> > dims = dimsizes(lat)
> > nlat = dims(0)
> > nlon = dims(1)
> >
> > gsn_define_colormap(wks,"WhBlGrYeRe")
> > res                         = True
> > res at gsnMaximize             = False
> > res at gsnDraw                 = False        ; Turn off for now.
> > res at gsnFrame                = False        ; Will draw later
> > res at tiMainOn                = False
> >
> > mpres                              = res
> > mpres at mpProjection                 = "LambertConformal"
> > mpres at mpLimitMode                  = "Corners"
> > mpres at mpLeftCornerLatF             = lat(0,0)
> > mpres at mpLeftCornerLonF             = lon(0,0)
> > mpres at mpRightCornerLatF            = lat(nlat-1,nlon-1)
> > mpres at mpRightCornerLonF            = lon(nlat-1,nlon-1)
> > mpres at mpLambertParallel1F          = b at TRUELAT1
> > mpres at mpLambertParallel2F          = b at TRUELAT2
> > mpres at mpLambertMeridianF           = b at STAND_LON
> > mpres at mpDataBaseVersion            = "MediumRes"
> > mpres at mpOutlineBoundarySets        = "GeophysicalAndUSStates"
> > mpres at mpFillOn                     = False
> > mpres at mpOutlineOn                  = True
> > mpres at mpCountyLineThicknessF       = 2.5
> > mpres at mpGeophysicalLineThicknessF  = 2.5
> > mpres at mpNationalLineThicknessF     = 2.5
> > mpres at mpProvincialLineThicknessF   = 2.5
> > mpres at mpUSStateLineThicknessF      = 2.5
> > mpres at mpCountyLineColor            = "Black"
> > mpres at mpGeophysicalLineColor       = "Black"
> > mpres at mpNationalLineColor          = "Black"
> > mpres at mpUSStateLineColor           = "Black"
> >
> > cnres                       = res
> > cnres at tfDoNDCOverlay        = False ; Testing
> > cnres at cnFillOn              = True
> > cnres at lbLabelBarOn          = True
> > cnres at cnLinesOn             = False
> > cnres at lbOrientation         = "Vertical"
> > cnres at sfDataArray           = ash
> > cnres at sfXArray              = lonN
> > cnres at sfYArray              = latN
> > cnres at sfMissingValueV       = ash at _FillValue
> > cnres at cnLevelSelectionMode  = "ManualLevels"
> > cnres at cnFillMode            = "RasterFill"
> > cnres at cnRasterSmoothingOn   = False
> > cnres at cnFillOpacityF        = 0.75
> > cnres at cnMinLevelValF        = 0.
> > cnres at cnMaxLevelValF        = 2.0
> > cnres at cnLevelSpacingF       = 0.05
> >
> > map         = gsn_csm_map(wks,mpres)
> > ash_contour = gsn_csm_contour(wks,ashNN,cnres)
> >
> > overlay(map,ash_contour)
> > draw(map)
> > 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/20160908/322bbbb7/attachment.html 


More information about the ncl-talk mailing list