[ncl-talk] how to plot observation stations wind vector on a map?

[dyjbean at gmail.com]

hi,
 i want add several obervation wind vector data on a map, but i cannot find appropriate function to plot.
 the following is my script:

++++++++++++++++++++++++++++++++++++++
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"
load "$GEODIAG_ROOT/geodiag.ncl"

begin
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; uv data
bj_csvs=systemfunc("ls uvdata/112920_inside_obs.csv")

lines=asciiread(bj_csvs,-1,"string")
nlines=dimsizes(lines) - 1

delim=","
field_names=str_split(lines(0),delim)
nfields=dimsizes(field_names)

;
fields=new((/nfields,nlines/),string)
do nf=0,nfields-1
  fields(nf,:)=str_get_field(lines(1:),nf+1,delim)
end do

print(fields(1,:) + "   "+fields(2,:) + "   "+fields(3,:)+"   "+fields(4,:))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; pm2.5 data

ncol = 3                                ; number of columns is 3£»ÄãµÄÁÐÊý
  stationdata= readAsciiTable("pm2.5_txt/15112920.txt", ncol, "float", 1)
  print(dimsizes(stationdata))
  nps  =dimsizes(stationdata)
  npts = nps(0)            ; Number of points.

  lat = stationdata(:,1)     ; latitude values
  lon = stationdata(:,0)     ; longitude values
  pm25= stationdata(:,2)   ; station numbers to appear on the map£»Õâ¸öµØ·½Äã¿ÉÒÔдÎÛȾÊý¾Ý

  printMinMax(pm25,0)
;-------interpolate-----------------------------------------------------
  olon =new(80,"float")
  olat =new(80,"float")
  data1=new((/80,80/),"float")
  do j=0,79
    olon(j)=113+j*0.1
  end do
  do k=0,79
    olat(k)=36+k*0.1
  end do
   


;---------set interpolate res--------------------------------------------
  olon!0        ="lon"
  olon at long_name="lontitude"
  olon at units    ="degrees-east"
  olon&lon      =olon

  olat!0        ="lat"
  olat at long_name="latitude"
  olat at units    ="degrees-north"
  olat&lat      =olat

;--------Call interpolation function--------------------------------------
  ;pm25 at _FillValue=999999.000000
  rscan=(/10,5,3,2,1,0.5,0.2,0.1/)
  data1=obj_anal_ic_deprecated(lon, lat, pm25, olon, olat, rscan, False);creanm
  printMinMax(data1 ,0)
  data1!0="lat"
  data1!1="lon"
  data1&lat=olat
  data1&lon=olon
  data1 at units=""
  data1 at long_name="pm2.5"
  ;data1 at _FillValue = pm25 at _FillValue
  printMinMax(olat,0)
  printMinMax(olon, 0)

;; plot
  type="png"
  wks = gsn_open_wks(type,"2015112920")     ; Open a workstation and
  colors=(/"White","Black","(/.027,.729,.145/)","(/.431,.824,.035/)",\
          "(/.800,.937,.012/)","(/1.00,.969,0.00/)","(/1.00,.871,0.00/)",\
          "(/1.00,.710,0.00/)","(/1.00,.549,0.00/)","(/0.00,.757,1.00/)",\
          "(/0.00,0.455,1.00/)","(/0.00,.122,1.00/)","(/1.00,.275,0.00/)",\
          "(/1.00,0.059,0.00/)"/)  
  gsn_define_colormap(wks,colors)

  res = True
  res at cnFillOn = True
  res at cnLinesOn = False
  res at cnLineLabelsOn = False
  res at gsnDraw                     = False
  res at gsnFrame                    = False
  res at gsnMaximize                 = True
  res at gsnAddCyclic = False

  
  res at mpLimitMode = "LatLon"
  res at mpMinLatF = 39.40
  res at mpMaxLatF = 41.00
  res at mpMinLonF = 115.35
  res at mpMaxLonF = 117.55

  res at isShowProvince =  False
  res at isShowBeijingCounty = True
  res at beijingCountyColor = "black"
  res at beijingCountyThickness = 1.0
  res at isShowBeijingRing2 = True
  res at beijingRing2Color = "black"
  res at beijingRing2Thickness = 1.50
  res at isShowBeijingRing3 = True
  res at beijingRing3Color = "black"
  res at beijingRing3Thickness = 1.50
  res at isShowBeijingRing4 = True
  res at beijingRing4Color = "black"
  res at beijingRing4Thickness = 1.5
  res at isShowBeijingRing5 = True
  res at beijingRing5Color = "black"
  res at beijingRing5Thickness = 1.5
  res at isShowBeijingRing6 = True
  res at beijingRing6Color = "black"
  res at beijingRing6Thickness = 1.5
  setup_china_map(res)

  plot = gsn_csm_contour_map(wks,data1,res)

  attach_china_map(wks, plot)

  draw(plot)

  ; wmsetp("vch", 0.05)
  ; wmsetp("vcc", 1)
  ; wmsetp("vcw", 3.)
  ; wmsetp("wbs", 0.05)
  ; wmsetp("col", 1)
  ; wmsetp("wbs", 0.06)
  ; wmsetp("ezf",2)
  ; wmsetp("wbs",0.12)
  wmsetp("ars",0.05)
  wmsetp("arl",1.5)
  ; wmvectmap(wks, tofloat(fields(2,:)), tofloat(fields(1,:)), tofloat(fields(3,:)), tofloat(fields(4,:)))
  ; wmvlbl(wks, 1.0, 0.0)
  wmbarbmap(wks, tofloat(fields(2,:)), tofloat(fields(1,:)), tofloat(fields(3,:)), tofloat(fields(4,:)))

  frame(wks)


end

++++++++++++++++++++++++++++++++++++++ 
the red and bold font is for wind barb.

the figure is as follows:
   
the wind barb are so small on this plot, but those sizes are so tiny , which function can be used to  enlarge the wind rod  for appropriate size.

thanks  



dyjbean at gmail.com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20161105/ede2379b/attachment-0001.html 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: image/png
Size: 77393 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20161105/ede2379b/attachment-0001.png