Masking data lying outside polyline

*I am using following program to plot contour using sfxarray . I do
not want to draw contour out side polyline .which function should i
used to masking data outside polyline.*
; Example script to produce plots for a WRF real-data run,
; with the ARW coordinate dynamics option.
; Script show how to zoom into a given area

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

begin

; Read the file
;************************************************

  lm=asciiread("lower_yamuna.prn",-1,"string")
  um=asciiread("upper_yamuna.prn",-1,"string")
  bt=asciiread("sahibi.prn",-1,"string")
; d1=asciiread("d1.prn",-1,"string")

 delim = " "
 nfields = str_fields_count(lm(1), delim)
 print(nfields)
 print(lm(1))

  lm_lat =stringtofloat(str_get_field(lm,nfields,delim))
  um_lat =stringtofloat(str_get_field(um,nfields,delim))
  bt_lat =stringtofloat(str_get_field(bt,nfields,delim))

  lm_lon =stringtofloat(str_get_field(lm,nfields-1,delim))
  um_lon =stringtofloat(str_get_field(um,nfields-1,delim))
  bt_lon =stringtofloat(str_get_field(bt,nfields-1,delim))

;***************************************************
;Read the input.csv file

  data=asciiread("input.csv",-1,"float")

  dim_lm=dimsizes(data)/3
 lon=new((/dim_lm/),"float")
 lat=new((/dim_lm/),"float")
 rh=new((/dim_lm/),"float")

 do i=0,dimsizes(data)-1,3
    lat(i/3)=data(i)
  end do
  do i=1,dimsizes(data)-1,3
    lon(i/3)=stringtofloat(data(i))
  end do

  do i=2,dimsizes(data)-1,3
    rh(i/3)=stringtofloat(data(i))
  end do

print(lat)
print(lon)
print(rh)

; For CSV file file following method is used

; delim = ","
 ;nfields = str_fields_count(lm(1), delim)
 ;print(nfields)
 ;print(lm(1))

  ;name =str_get_field(lm,nfields-3,delim)

 ; lat =stringtofloat(str_get_field(lm,nfields-2,delim))

 ; lon =stringtofloat(str_get_field(lm,nfields-1,delim))

 ; rh = stringtofloat(str_get_field(lm,nfields,delim))

k=0

do i=0,dimsizes(lat)-1

if (.not.any(ismissing(rh(i)))) then

 rh(k)=rh(i)

 lon(k)=lon(i)

 lat(k)=lat(i)

 k=k+1

 end if

 end do
;******************************************************
;save station name,lat ,lon,rainfall values into an array.

rhn=new(k,float)
lond=new(k,float)
latd=new(k,float)

do i=0,k-1

rhn(i)=stringtofloat(rh(i))
lond(i)=stringtofloat(lon(i))
latd(i)=stringtofloat(lat(i))
 end do

print(rhn)
print(lond)
print(latd)
;********************************************************

;setting some plotting option & contour plot.

  wks = gsn_open_wks("pdf","contour")

  colors = (/"white","black",
"white","yellow","yellowgreen","forestgreen","deepskyblue","royalblue4","darkorange3","red"
/)

  gsn_define_colormap(wks,colors)

   res = True
   res@gsnMaximize = True
   res@gsnDraw=False
   res@gsnFrame=False
   res@cnLevelSelectionMode = "ExplicitLevels" ; set explicit contour levels
   res@cnLevels = (/1.,10.,20.,40.,70.,130.,200./)
   res@cnFillOn =True
   res@cnLinesOn =False;True
   res@cnFillMode ="AreaFill"
   res@cnSmoothingOn =True
; res@cnSmoothingDistanceF =.0000001
 ; res@cnSmoothingTensionF =5.5
  ; res@cnMaxPointDistanceF =0
 ; res@cnFillScaleF = 9

   res@mpProjection = "Mercator" ; Change the map projection.
  ; res@cnRasterSmoothingOn =True
  res@lbLabelAutoStride = True
  res@lbBoxLinesOn = False
  res@tiMainString = "Mahanadi Basin (Rainfall in mm)"
 res@gsnAddCyclic=True; False
; res@tfDoNDCOverlay=True
  res@sfXArray = lond
  res@sfYArray = latd

   res@mpLimitMode = "LatLon" ; Change the area of the map
    res@mpMinLatF = 27.0 ; viewed.
    res@mpMaxLatF = 31.5
    res@mpMinLonF = 75.5
    res@mpMaxLonF = 78.8

 ; res@mpLimitMode = "Corners"
  ;res@mpLeftCornerLatF =23;min(lm_lat)-.1;min(latd)-.1;lat2d(20,18)
  ;res@mpLeftCornerLonF =72;min(um_lon)-.1;min(lond)-1.9
;lon2d(20,18)
  ;res@mpRightCornerLatF =35;max(um_lat)+.1;max(latd)+.3; lat2d(58,52)
  ;res@mpRightCornerLonF =90;max(lm_lon)+.1;max(lond)+.3;lon2d(58,52)

 res@mpOutlineDrawOrder = "PostDraw"
  res@mpFillDrawOrder = "PreDraw"
  res@pmTickMarkDisplayMode = "Always"
   plot = gsn_csm_contour_map(wks,rhn,res)
; plot = gsn_csm_contour_map_ce(wks,rhn,res)

;*********************************************************************
; Now add some markers to show where the original 1D points are.
;
  ;mkres = True
  ;mkres@gsMarkerIndex = 16 ; Filled dots
  ; mkres@gsMarkerSizeF = .001
 ;dum = gsn_add_polymarker(wks,map,lond,latd,mkres)

;*******************************************************
;draw basin bounadry

 ; draw basin bounadry

  poly_res=True
 ; poly_res@gsLineColor="Red"
  poly_res@gsLineThicknessF=2.5

plot1=gsn_add_polyline(wks,plot,lm_lon,lm_lat,poly_res)
plot2=gsn_add_polyline(wks,plot,um_lon,um_lat,poly_res)
plot3=gsn_add_polyline(wks,plot,bt_lon,bt_lat,poly_res)

  txres = True
  txres@txFontHeightF = 15;fheight ; Set the font height
; txres@txConstantSpacingF=1.2

; gsn_text_ndc(wks,"Sub basin wise Average Rainfall Forecast",.4,.27,txres)

  delete(txres)
  txres = True
  txres@txFontHeightF = 12;fheight ; Set the font height
; txres@txConstantSpacingF=1.2

   plot9=gsn_add_text(wks,plot,"Lower Yamuna",77.2,29,txres)
  plot9=gsn_add_text(wks,plot,"Upper Yamuna",77.7,31,txres)
  plot9=gsn_add_text(wks,plot,"Sahibi",76.7,28,txres)

  draw(plot)
  frame(wks)

end

-- 
Abhishek Kumar Gupta,
Junior Research Fellow,
India Meteorological Department,
Lodi Road ,New Delhi-110003,
India

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Tue Mar 6 23:37:24 2012