RE: [ncl-talk] random extraction of data

From: Dennis Shea <shea_at_nyahnyahspammersnyahnyah>
Date: Fri, 2 Mar 2007 10:52:48 -0700 (MST)

>Thanks for your solution. It is working.
>But iu takes the maximum value 3200 whereas I have 10500 grid points.
>I did not figure it out how will I increase the value of iu (it will be good
>if iu takes the value nearer to 9000) so that random 500 grids will be more
>disperse over the region (i.e. 10500 grid points)
>
>-----Original Message-----

I took another look. Attached is a simple function.
This works fine for small grids. Otherwise, it may be
a bit slow.

Please note: Just because the grid locations are randomly
             sampled does not mean that the sampled values
             are random (independent).
             
             ie: if (say) the 37th random subscript
                 corresponds to T(15,38) and the (say)
                 437th random subscript corresponds
                 to T(15,39), clearly T(15,38)
                 and T(15,39) are not independent.
             

       

function unique_index (nMax[1], nUnique[1]:integer, nGenerate[1]:integer)
local nRandom, ratio, IR, LU, n, i, ni, iu, niu
begin
   if (nUnique.gt.nMax) then
       print("-------")
       print("unique_subscripts: nUnique > nMax ; impossible")
       print("-------")
       return(-999)
   end if

; to ensure that enough random integers are generated check the ratio
     
   nRandom = nMax
   ratio = nRandom/nUnique
   if (ratio.lt.5 .or. nGenerate.gt.5) then
       nRandom = nRandom*max( (/5,nGenerate/) )
   end if

   IR = round(random_uniform(0,(nMax-1),nRandom), 3)

   LU = new( nRandom, "logical")
   LU = True
                     ; this *will* be slow if nRandom is large
   do n=0,nRandom-1
      i = ind(IR.eq.IR(n) .and. LU(n))
      ni = dimsizes(i)
      if (ni.gt.1) then
         LU(i(1:)) = False
      end if
      delete(i)
   end do
 
   iu = ind(LU) ; index of unique subscripts
   niu = dimsizes(iu) ; number of unique subscripts
   if (niu.ge.nUnique) then
       return(IR(iu(0:nUnique-1))) ;
   else
       print("-------")
       print("unique_subscripts: no. of unique subscripts < NR ; increase nGenerate")
       print("-------")
       return(IR(iu(0:niu-1)))
   end if
       
end

begin
   nlat = 70
   mlon = 150
   nUnq = 500 ; number of unique subscripts

                               ; generate a 2D array
   T = generate_2d_array(10,10,-10,40,0,(/nlat,mlon/))

                               ; product(dimsizes(T)) = nlat*mlon
   iT = unique_index (product(dimsizes(T)), nUnq, 0)

   T1D = ndtooned( T )
   
   T_rdmspl = T1D( iT ) ; randomly sampled T
   print(iT+" "+T_rdmspl)
   
; or
   print("=======================")
   ir = ind_resolve(iT, dimsizes(T) ) ; 2D [500,2]
   printVarSummary(ir)

   do n=0,nUnq-1
      print("T: "+T(ir(n,0),ir(n,1))+" ir: "+ir(n,0)+" "+ir(n,1))
   end do
end

_______________________________________________
ncl-talk mailing list
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Fri Mar 02 2007 - 10:52:48 MST

This archive was generated by hypermail 2.2.0 : Fri Mar 02 2007 - 17:25:16 MST