I did this Monday evening and, inadvertently, did not post
==========================================================
If 'a' and 'b' are one dimensional
c = a ; same dimensionality as 'a'
c = 0.0 ; array syntax (init to 0.0)
ii = ind(a.gt.0.0)
c(ii) = b(ii)/a(ii)
---- A slight variation on Kyle's suggestion is the following. Note: This approach has a bug which has been fixed in NCL v6.2.0 so it can not be used thru 6.1.2 c = 1.0/where(a .gt. 0, b/a, a@_FillValue) This assumes 'a' has an _FillValue 'c' will have _FillValue in all instances where a.le.0 To get 0s, you would have to do the following on the next lines c@_FillValue = 0.0 ; all _FillValue set to 0.0 delete(c@_FillValue) A possible unfortunate side effect would be cases where 'a' had_FillValue prior to the calculation. Then it would be impossible to distinguish between the original suite of a@_FillValue and the above values. --- quick and dirty 'solution' undef("my_zero_where") function my_zero_where(a,b) ; a can be any dimensionality ; rank[a] = rank[b] local dima, dimb, ranka, rankb, c, a1d, b1d, c1d, ii, jj, iir, jjr, n\ , dimiir, dimjjr, niir, njjr begin dima = dimsizes(a) ranka = dimsizes(dima) dimb = dimsizes(b) rankb = dimsizes(dimb) if (ranka.ne.rankb) then print("my_zero_where: a,b do not conform") exit end if if (all(a.gt.0.0)) then ; nothing special need be done c = b/a ; array syntax ... 'a' any dimensionality else c = (/ a /) ; create array; same a@_FillValue if (ranka.eq.1) then ii = ind(a.gt.0.0) c(ii) = b(ii)/a(ii) jj = ind(a.le.0.0) c(jj) = 0.0 else a1d = ndtooned(a) b1d = ndtooned(b) c1d = (/ a1d /) ii = ind(a1d.gt.0.0) c1d(ii) = b1d(ii)/a1d(ii) jj = ind(a1d.le.0.0) c1d(jj) = 0.0 ; 'resolve' [remap] 1d indices to to nD iir = ind_resolve(ii,dima) dimiir= dimsizes(iir) niir = dimiir(0) ; number of elements a1d > 0.0 if (ranka.ne.3) then print("my_zero_where: the nD part supports 3D only") print(" The current array is "+ranka+"D") exit end if do n=0,niir-1 c(iir(n,0),iir(n,1),iir(n,2)) = c1d(ii(n)) end do jjr = ind_resolve(jj,dima) dimjjr= dimsizes(jjr) ; number of elements a1d <= 0.0 njjr = dimjjr(0) ; number of elements a1d > 0.0 do n=0,njjr-1 c(jjr(n,0),jjr(n,1),jjr(n,2)) = c1d(jj(n)) end do end if end if return(c) end ;++++++++++++++++++++ TEST DRIVER ++++++++++++++++++++ ; MAIN N = 100 A1D = random_normal(0,5,N) B1D = random_uniform( -5, 5, N) A1D@_FillValue = 1e20 IX = ispan(0,N-1,5) A1D(IX) = A1D@_FillValue I0 = ispan(3,N-1,5) A1D(I0) = 0.0 C1D = my_zero_where(A1D,B1D) print(C1D) On 1/20/14, 7:06 AM, Robert Schuster wrote: > Thank you! Working with missing values is a good idea! > It is good to know, that the expressions for both cases are evaluated for all values in the arrays and not only for those where the condition is true. > > Robert > > Am 20.01.2014 um 14:41 schrieb Kyle Griffin <ksgriffin2@wisc.edu>: > >> Unfortunately, Robert, that is exactly how the where function is defined to work. I suggest turning the zeros into missing values (at least temporarily) for the calculation and then subbing the values from the original function back in with another where statement, but others may have better suggestions to work around this issue. >> >> >> Kyle >> >> ---------------------------------------- >> Kyle S. Griffin >> Department of Atmospheric and Oceanic Sciences >> University of Wisconsin - Madison >> Room 1421 >> 1225 W Dayton St, Madison, WI 53706 >> Email: ksgriffin2@wisc.edu >> >> >> On Mon, Jan 20, 2014 at 3:05 AM, Robert Schuster <rxschuster@gmail.com> wrote: >> Hello, >> >> I have two arrays that I want to divide by each other, one of them contains zeros: >> >> a = (/0,1/) >> b = (/1,1/) >> >> I tried to use the where function to avoid a divide by zero error: >> >> c = where(a .gt. 0, b/a, 0) >> >> But interestingly the divide operation is performed on the whole arrays and not only where the condition a .gt. 0 is true: >> >> fatal:divide: Division by 0, Can't continue >> fatal:Div: operator failed, can't continue >> >> Is that the way the where function is supposed to work? >> >> Robert >> _______________________________________________ >> ncl-talk mailing list >> List instructions, subscriber options, unsubscribe: >> http://mailman.ucar.edu/mailman/listinfo/ncl-talk >> > > > > > _______________________________________________ > ncl-talk mailing list > List instructions, subscriber options, unsubscribe: > http://mailman.ucar.edu/mailman/listinfo/ncl-talk > _______________________________________________ ncl-talk mailing list List instructions, subscriber options, unsubscribe: http://mailman.ucar.edu/mailman/listinfo/ncl-talkReceived on Wed Jan 22 10:44:15 2014
This archive was generated by hypermail 2.1.8 : Fri Feb 07 2014 - 16:39:11 MST