[ncl-talk] attaching a shapefile to a wrfoutput plot

[Kerandi, Noah Misati (IMK)]

Dear all,

I have tried one example of ncl script for plotting surface wind (see below) to attach a shape file to  a wrf output plot however with no success.

The shape file doesn't get  attached to the two plots.

I will be glad if somebody moves me out from this state.

Thank you.


Noah
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
;load "./WRFUserARW.ncl"

begin
; ; This needs to have a ".nc" appended, so just do it.
  a = addfile("file","r")



; We generate plots, but what kind do we prefer?
   type = "x11"
  ;type = "pdf"
  ;type = "ps"
; type = "ncgm"
  wks = gsn_open_wks(type,"~/test")

; Set some basic resources
  res = True

  res at MainTitle                   = "Surface wind 10m"


  mpres = True

 ; Set some Basic Plot options
  res = True


  res at Footer = False

  mpres = True
 pltres            = True   ; Basic overlay plot options


  mpres at mpNationalLineColor         = "Black"
  mpres at mpUSStateLineColor          = "Black"
  ;mpres at mpGridLineColor             = "Black"
  mpres at mpLimbLineColor             = "Black"
  mpres at mpPerimLineColor            = "Black"
  mpres at mpGeophysicalLineThicknessF = 2.5
  ;mpres at mpGridLineThicknessF       = 2.0
  ;mpres at mpLimbLineThicknessF       = 2.5
  mpres at mpNationalLineThicknessF    = 2.5
   mpres at mpDataBaseVersion           = "MediumRes"

  mpres at mpFillBoundarySets          ="AllBoundaries"
  mpres at mpOutlineBoundarySets       = "National"

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; What times and how many time steps are in the data set?
  times = wrf_user_getvar(a,"times",-1)  ; get all times in the file
  ntimes = dimsizes(times)         ; number of times in the file
shp = new(ntimes,graphic)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  do it = 0,ntimes-1,2             ; TIME LOOP

    print("Working on time: " + times(it) )
    res at TimeLabel = times(it)   ; Set Valid time to use on plots


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; First get the variables we will need

    slp = wrf_user_getvar(a,"slp",it)    ; slp
      wrf_smooth_2d( slp, 3 )            ; smooth slp
    tc = wrf_user_getvar(a,"tc",it)      ; 3D tc
    td = wrf_user_getvar(a,"td",it)      ; 3D td
    u  = wrf_user_getvar(a,"ua",it)      ; 3D U at mass points
    v  = wrf_user_getvar(a,"va",it)      ; 3D V at mass points
    td2 =  wrf_user_getvar(a,"td2",it)   ; Td2 in C
    tc2 = wrf_user_getvar(a,"T2",it)     ; T2 in Kelvin
       tc2 = tc2-273.16                  ; T2 in C
    u10 = wrf_user_getvar(a,"U10",it)    ; u at 10 m, mass point
    v10 = wrf_user_getvar(a,"V10",it)    ; v at 10 m, mass point
     rh = wrf_user_getvar(a,"rh",it)        ;
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


     ;tf2 = 1.8*tc2+32.                    ; Turn temperature into Fahrenheit
      tc2 at description = "Surface Temperature"
      tc2 at units = "~F34~0~F~C"
    ;td_f = 1.8*td2+32.                   ; Turn temperature into Fahrenheit
      td2 at description = "Surface Dew Point Temp"
      td2 at units = "~F34~0~F~C"
    ;u10 = u10*1.94386                    ; Turn wind into knots
    ;v10 = v10*1.94386
      u10 at units = "m/s"
      v10 at units = "m/s"


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

    ; Plotting options for T
      opts = res
      opts at cnFillOn = True
      opts at ContourParameters = (/ 0., 30., 5./)
      opts at gsnSpreadColorEnd = -3  ; End third from the last color in color map
      contour_tc = wrf_contour(a,wks,tc2,opts)

      delete(opts)


    ; Plotting options for Td
      opts = res
      opts at cnFillOn = True
      opts at cnLinesOn = True
      opts at cnLineLabelsOn = True
      opts at ContourParameters = (/ 0., 30., 5./)
      opts at cnLineLabelBackgroundColor = -1
      opts at gsnSpreadColorEnd = -3  ; End third from the last color in color map
      contour_td = wrf_contour(a,wks,td2,opts)

      delete(opts)


    ; Plotting options for SLP
      opts = res
      ;opts at cnLineColor = "Blue"
      ;opts at cnHighLabelsOn = True
      ;opts at cnLowLabelsOn = True
       opts at ContourParameters = (/ 1000., 1200., 2. /)
      ;opts at cnLineLabelBackgroundColor = -1
      ;opts at gsnContourLineThicknessesScale = 2.0
      contour_psl = wrf_contour(a,wks,slp,opts)

      delete(opts)


   ; Plotting options for Wind Vectors

    opts = res
        opts at FieldTitle = "Wind"   ; overwrite Field Title
        opts at NumVectors = 47       ; wind barb density
        opts at vcRefAnnoOrthogonalPosF = 0.02  ; move ref vector up
        opts at vcWindBarbColor  = "black"       ; Draw wins barbs in white
        opts at vcRefMagnitudeF  = 10.0          ; define vector ref mag
        opts at vcRefLengthF     = 0.045            ; define length of vec ref
        opts at vcGlyphStyle     = "CurlyVector"    ; turn on curly vectors
        opts at vcMinDistanceF    = 0.055      ; thin out windbarbs

        opts at vcRefAnnoOn       =True
        opts at gsnFrame          = False            ; so we can draw time stamp
        opts at vcRefAnnoSide     = "Right"
         vector = wrf_vector(a,wks,u10,v10,opts)

         delete(opts)


    ; MAKE PLOTS

     plot = wrf_map_overlays(a,wks,(/contour_tc,contour_psl,vector/),pltres,mpres)


     plot= wrf_map_overlays(a,wks,(/contour_td,vector/),pltres,mpres)
 ;define shapefiles
shp_filename = ("./shape.shp")
;-- set shapefile resources
  shpres                    =  True
  shpres at gsLineThicknessF   =  3.9                   ;-- increase line thickness
  shpres at gsLineColor        = "Purple"             ;-- line colorgsLineThicknessF

shp= gsn_add_shapefile_polylines(wks,plot,shp_filename,shpres)
 draw(plot)
frame(wks)
end do ; END OF TIME LOOP

end

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20160414/cd726797/attachment.html