[ncl-talk] Overlaying Contour on WRF Domain

Genito Amos Maure genito.maure at uem.mz
Thu Sep 8 09:57:37 MDT 2016


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





More information about the ncl-talk mailing list