Example 1 - find nearest points to a given point in 3-space (NCL)


load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_code.ncl"

begin
;
; Initialize some work arrays.
;
  mtri = 150000
  rtri = new((/mtri,10/),float)
  rtwk = new(2*mtri,float)
  itwk = new(mtri,integer)

;
; Create the input arrays that will hold the indices of the points we
; want to plot.
;
  n        = 1331
  nearest  = 500
  nfarther = n - nearest

  npts = new(nearest,integer)
  mpts = new(nfarther,integer)
  
;
; Generate arrays of randomly-positioned points in the unit cube.
;
  x = random_uniform(0.,1.,n)
  y = random_uniform(0.,1.,n)
  z = random_uniform(0.,1.,n)

;
; Specify the reference point from which we want to find the nearest
; "nearest" points.
;
  px = 0.5
  py = 0.5
  pz = 0.5

;
; Open a workstation.
;
  wks = gsn_open_wks("x11","scatter")
;
; Initialize Tdpack parameters.
;
; Move the plot up a bit.
;
  tdsetp("VPB",0.09)
  tdsetp("VPT",0.99)
  tdsetp("VPL",0.11)
  tdsetp("VPR",1.00)
  tdinit ((/4.6,  3.0, 3.3/), (/0.5, 0.5, 0.5/), (/0.5, 0.5, 2.7/),0.)
;
; Set up some colors.
;
  tdclrs (wks, 1, 0., 0.8, 8, 37, 8)
;
; Draw some labels so we can see them while waiting for plot to draw.
;
  txres = True
  txres@txFontHeightF = 0.02
  gsn_text_ndc(wks,":F22:Find the nearest N points in three space", \
                   0.50,0.95,txres)

  txres@txJust = "CenterLeft"
  gsn_text_ndc(wks,":F21:Red ball = reference point", \
                   0.01,0.17,txres)
  gsn_text_ndc(wks,":F21:Green balls = near points", \
                   0.01,0.12,txres)
  gsn_text_ndc(wks,":F21:Gray balls = far points", \
                   0.01,0.07,txres)
;
; Define style indices for shades of gray, green, and red.
;
  tdstrs (1,  8, 37,   8,  37, 1, 1, 0, 0.05, 0.05, 0.)
  tdstrs (3,  8, 37,  68,  97, 1, 1, 0, 0.05, 0.05, 0.)
  tdstrs (4,  8, 37,  98, 127, 1, 1, 0, 0.05, 0.05, 0.)
;
; Store the indices of the nearest points in NPTS and the complement
; of that set (with respect to the entire input dataset) in MPTS.
;
  npts(0) = shgetnp(px,py,pz,x,y,z,0)
  do i=2,nearest
    npts(i-1) = shgetnp(px,py,pz,x,y,z,1)
  end do
  do i=0,n-nearest-1
    mpts(i) = shgetnp(px,py,pz,x,y,z,1)
  end do
;
; Add near points in green.
;
  ntri = 0
  dotsize = 0.02
  do i=0,nearest-1
    inx = npts(i)
    tdmtri(-5,x(inx),y(inx),z(inx),dotsize,rtri,ntri,4,0.,0.,0.,1.,1.,1.)
  end do
;
; Add the farther points in gray.
;
  do i=0,nfarther-1
    inx = mpts(i)
    tdmtri(-5,x(inx),y(inx),z(inx),dotsize,rtri,ntri,1,0.,0.,0.,1.,1.,1.)
  end do
;
; Mark the reference point in red.
;
  tdmtri(-5,px,py,pz,1.2*dotsize,rtri,ntri,3,0.,0.,0.,1.,1.,1.)
;
;  Draw dots.
;
  tdotri(rtri,ntri,rtwk,itwk,0) 
  tddtri(wks,rtri,ntri,itwk) 
;
; Draw a box around the perimeter.
;
  tdgrds (wks, (/0., 1., 0./), (/1., 0., 1./), (/-1., -1., -1./),11,0)
  tdgrds (wks, (/0., 1., 0./), (/1., 0., 1./), (/-1., -1., -1./),11,1)

  frame(wks)
end

home | contents | defs | procedures | examples | errors