# 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