Lat Lon plot from WRF grib file

Hi,

I want to make a plot from a WRF grib file (obtained from WPP).

I'm trying to change the coordinate system from x/y to lat/lon, but I get this error message:

warning:tmXBMinorLengthF is not a valid resource in 0_contour.PlotManager at this time
warning:NhlGetValues:Error retrieving tmXBMinorLengthF

I need to overlay a shapefile which is in lat lon and that is why I have to plot the WRF output in lat lon.

I think the problem is that I'm not doing the conversion properly.

Thanks for your help.

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/csm/contributed.ncl"
;************************************************
; MAIN
;************************************************
begin

DATADir = "/home/estatio/ncl_scripts/pbl/"
  FILES = systemfunc (" ls -1 " + DATADir + "WRF* ")
  numFILES = dimsizes(FILES)

do ifil = 0,numFILES-1

a=addfile(FILES(ifil)+".grb","r")

pbl=a->HPBL_GDS3_SFC_10
lat2d=a->NLAT_GDS3_SFC_10(:,:)
lon2d=a->ELON_GDS3_SFC_10(:,:)
pltType = "x11"

 wks = gsn_open_wks (pltType,ifil) ; open workstation
    gsn_define_colormap(wks,"BlueDarkRed18") ; choose colormap

    res = True ; plot mods desired
    res@gsnMaximize = True ; make ps, eps, pdf large
    res@gsnSpreadColors = True ; Use full color map
 ; res@gsnSpreadColorEnd= -2 ; do not use gray for contours
    res@cnFillOn = True ; color
    res@cnLinesOn = False
    res@cnLineLabelsOn = False

    res@lbLabelAutoStride= True ; let NCL choose stride

    res@gsnDraw = False ; don't draw
    res@gsnFrame = False

 res@tiXAxisString = lon@long_name
 res@tiYAxisString = lat@long_name
 res@sfXArray = lon2d
 res@sfYArray = lat2d
 
plot = gsn_csm_contour(wks,pbl,res) ; Pacific

; Open the states shapefile and draw the state boundaries...
  f = addfile("/home/estatio/ncl_scripts/nybb_09c_av/nybb.shp", "r")

  geometry = f->geometry
  segments = f->segments
  segsDims = dimsizes(segments)
  dims = dimsizes(geometry)

;
; We need to reference these inside a nested do loop, so it's
; better to read them off the file here once, and use the
; local variables inside the do loops.
;
  geom_segIndex = f@geom_segIndex
  geom_numSegs = f@geom_numSegs
  segs_xyzIndex = f@segs_xyzIndex
  segs_numPnts = f@segs_numPnts

  plres = True
  plres@gsLineColor = "black"
  plres@gsLineThicknessF=1.5
  lines = new(segsDims(0),graphic)

  numFeatures = dims(0)

  segNum = 0
  latc = f->y
  lonc = f->x
  do i=0, numFeatures-1
     startSegment = geometry(i, geom_segIndex)
     numSegments = geometry(i, geom_numSegs)
     do seg=startSegment, startSegment+numSegments-1
        startPT = segments(seg, segs_xyzIndex)
        endPT = startPT + segments(seg, segs_numPnts) - 1
        lines(segNum) = gsn_add_polyline(wks, plot, lonc(startPT:endPT), \
                                                    latc(startPT:endPT), plres)
        segNum = segNum + 1
     end do
end do

draw(plot)
frame(wks)

 delete(wks)

; cmd = "convert -geometry 1200x1200 -density 1200 -trim " + ifil + ".ps " + \
; ifil + ".png"
; system(cmd)

end do

end

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Thu Mar 11 19:09:57 2010