[IDV #VVW-585732]: IDV - limits of visad?

Hi Tyn-

Sorry for the delay in responding.

> Institution: ITC
> Package Version: 2.3 build date:2007-08-15 19:04 UTC
> Operating System: Linux
> Hardware Information: Java: home: /usr/lib/jvm/ia32-java-6-sun-1.6.0.00/jre 
> version: 1.6.0 j3d:1.3.1
> Inquiry: Hi Jeff/Don,
> 
> i think i have approached the limits of visad, but hope i haven't ;-)
> 
> In the jython script below i'm looping over all coordinates (crd) and all 
> times (times),
> but the function:
> 
> org.apache.commons.math.analysis import SplineInterpolator, 
> PolynomialSplineFunction, PolynomialFunction
> 
> demands i can update a coordinates of a singelBandedImage (time)
> in a imageSequence without explicitely addressing the singleBandedImage.

I'm not sure what you are trying to do here exactly.  Do you want
to set the value on a particular pixel, or calculate all the values
and set the pixels.  

> So analog to:
> 
> currTimeSequence.getSample(time).getFloats(0)
> 
> i would like to be able to:
> 
> currTimeSequence.getSample(time).setFloats(0)

If you call getFloat(0), you should get the values without
copying them.  So, if you then modify the values in the
resulting arrays, it should modify the original image
values.

> without closing the loop-over-coordinates loop.

I'm still not clear what the above statement would do.  There is a 
setSamples(newvals,0)
method you could call to set new values in the function, but I'm not sure that's
waht you want.

> This would maintain the "status" of the function per pixel (which is a binary 
> object),
> as i cannot save a binary object in the pixel i have to keep it in memory.
> 
> Sorry for being so unclear, i'm not a real programmer....

I just play one on TV. ;-)

Don

> ##########################################################################################################################################################
> # compute the Cloud Index of timesequence imagery (call this jython function 
> only for image sequences)
> ##########################################################################################################################################################
> # Purpose: compute the Cloud Index of timesequence imagery (call this jython 
> function only for image sequences)
> # Creator: V. Venus
> #
> 
> 
> 
> 
> def 
> getCloudIndexHeliosatImageSeqN52(singleBand,histTimeSequence,histCsaSeq,currTimeSequence,currCsaSeq,user_pGround,user_pCloud):
> import sys;
> sys.add_package('visad.meteorology');
> sys.add_package('GridUtil');
> 
> sys.add_package('org.apache.commons.math.stat.descriptive.rank');
> from org.apache.commons.math.stat.descriptive.rank import Percentile
> 
> sys.add_package('org.apache.commons.math.analysis');
> from org.apache.commons.math.analysis import UnivariateRealInterpolator, 
> UnivariateRealFunction
> from org.apache.commons.math.analysis import SplineInterpolator, 
> PolynomialSplineFunction, PolynomialFunction
> 
> from java.util import TimeZone
> from jarray import zeros, array
> from java.lang import Float, Double, Object
> 
> from visad import DateTime
> from visad.meteorology import ImageSequence
> from visad.meteorology import SingleBandedImage
> 
> domain = singleBand.getDomainSet()
> cs = domain.getCoordinateSystem()
> len =  domain.getLength()
> 
> # Clone the incoming objects
> news = singleBand.clone()
> newd = currTimeSequence.clone()
> 
> alg = Percentile()
> pG = Double(user_pGround) # Sets the value of the quantile field (which 
> determines what percentile is computed when
> # evaluate() is called with quantile argument).
> 
> pC = Double(user_pCloud) # Sets the value of the quantile field (which 
> determines what percentile is computed when
> # evaluate() is called with quantile argument).
> 
> knot4 = zeros(13,"d")
> bin4 = zeros(13,"d")
> bin = zeros(13, Object)
> count = zeros(1,"d")
> 
> # get length of imageSequence
> seq = histTimeSequence.getDomainSet().getLength()
> currSeq = currTimeSequence.getDomainSet().getLength()
> # create array to store values
> for i in range(bin.__len__()):
> bin[i] = zeros(seq, "d")
> 
> for crd in xrange(len):
> count = zeros(bin.__len__(), "i")
> # get all temperatures for one coordinate for all times time and put them 
> into arr
> #print "Entering nested forloop 1"
> for time in xrange(seq):
> refl = histTimeSequence.getSample(time).getFloats(0)
> #print "    Value for current pixel (#", crd + 1, "): ",values[0][crd]
> CA = histCsaSeq.getSample(time).getFloats(0)
> refl = refl[0][crd]/100
> if (refl > 0.0075):
> if (CA[0][crd] < 10):
> #print "CA < 10"
> bin[0][count[0]] = refl
> #knot[0] = 0
> count[0] = count[0] + 1
> #print "CA: ",CA," refl: ", refl
> elif (CA[0][crd] >= 10) and (CA[0][crd] < 20):
> #print "(CA >= 10) and (CA < 20)"
> bin[1][count[1]] = refl
> #knot[1] = 10
> count[1] = count[1] + 1
> #print "CA: ",CA," refl: ", refl
> elif (CA[0][crd] >= 20) and (CA[0][crd] < 30):
> #print "(CA >= 20) and (CA < 30)"
> bin[2][count[2]] = refl
> #knot[2] = 20
> count[2] = count[2] + 1
> #print "CA: ",CA," refl: ", refl
> elif (CA[0][crd] >= 30) and (CA[0][crd] < 40):
> #print "(CA >= 30) and (CA < 40)"
> bin[3][count[3]] = refl
> #knot[3] = 30
> count[3] = count[3] + 1
> #print "CA: ",CA," refl: ", refl[0][i]
> elif (CA[0][crd] >= 40) and (CA[0][crd] < 50):
> #print "(CA >= 40) and (CA < 50)"
> bin[4][count[4]] = refl
> #knot[4] = 40
> count[4] = count[4] + 1
> #print "CA: ",CA," refl: ", refl
> elif (CA[0][crd] >= 50) and (CA[0][crd] < 60):
> #print "(CA >= 50) and (CA < 60)"
> bin[5][count[5]] = refl
> #knot[5] = 50
> count[5] = count[5] + 1
> #print "CA: ",CA," refl: ", refl
> elif (CA[0][crd] >= 60) and (CA[0][crd] < 70):
> #print "(CA >= 60) and (CA < 70)"
> bin[6][count[6]] = refl
> #knot[6] = 60
> count[6] = count[6] + 1
> #print "CA: ",CA," refl: ", refl
> elif (CA[0][crd] >= 70) and (CA[0][crd] < 80):
> #print "CA > 70"
> bin[7][count[7]] = refl
> #knot[7] = 70
> count[7] = count[7] + 1
> #print "CA: ",CA," refl: ", refl
> elif (CA[0][crd] >= 80) and (CA[0][crd] < 90):
> #print "CA > 80"
> bin[8][count[8]] = refl
> #knot[8] = 80
> count[8] = count[8] + 1
> #print "CA: ",CA," refl: ", refl
> elif (CA[0][crd] >= 90) and (CA[0][crd] < 100):
> #print "CA > 90"
> bin[9][count[9]] = refl
> #knot[9] = 90
> count[9] = count[9] + 1
> #print "CA: ",CA," refl: ", refl
> elif (CA[0][crd] >= 100) and (CA[0][crd] < 110):
> #print "CA > 100"
> bin[10][count[10]] = refl
> #knot[10] = 100
> count[10] = count[10] + 1
> #print "CA: ",CA," refl: ", refl
> elif (CA[0][crd] >= 110) and (CA[0][crd] < 120):
> #print "CA > 110"
> bin[11][count[11]] = refl
> #knot[11] = 110
> count[11] = count[11] + 1
> #print "CA: ",CA," refl: ", refl
> elif (CA[0][crd] >= 120) and (CA[0][crd] < 130):
> #print "CA > 120"
> bin[12][count[12]] = refl
> #knot[12] = 120
> count[12] = count[12] + 1
> #print "CA: ",CA," refl: ", refl
> knot = 0
> valid = 0
> binLength = bin.__len__()
> print "binlength: ", binLength
> for i in xrange(binLength):
> if (count[i] > 0):
> print "valid"
> binTemp = zeros(count[i], "d")
> for j in xrange (0,count[i]):
> binTemp[j] = bin[i][j]
> print "binTemp: ", binTemp
> bin4[valid] = alg.evaluate(binTemp, pG)
> knot4[valid] = knot
> valid = valid + 1
> else:
> print "empty bin at ", i
> #invalid = invalid + 1
> knot = knot + 10
> cnt = 0
> for i in xrange(5):
> if (count[i] > 0):
> cnt = cnt + count[i]
> binCloud = zeros(cnt, "d")
> #x = 0
> j = 0
> i = 0
> index = 0
> for i in xrange(5):
> if (count[i] > 0):
> for j in xrange (0,count[i]):
> binCloud[index] = bin[i][j]
> index = index + 1
> else:
> print "empty bin at ", i
> print "binCloud: ", binCloud
> rc = alg.evaluate(binCloud, pC)
> print "rc: ", rc
> 
> bin4Temp = zeros(valid, "d")
> knot4Temp = zeros(valid, "d")
> for j in xrange (0,valid):
> bin4Temp[j] = bin4[j]
> knot4Temp[j] =  knot4[j]
> print "knot4: ", knot4Temp
> print "bin4: ", bin4Temp
> knotLength = knot4Temp.__len__()
> print "knotLength: ", knotLength
> if (knotLength > 2):
> print "Sorry, at least 3 datapoints are required to compute a spline 
> interpolant. Try again using a longer timeseries."
> 
> interpolator = SplineInterpolator();
> function = interpolator.interpolate(knot4Temp,bin4Temp);
> for time in xrange(currSeq):
> valuesCSA = currCsaSeq.getSample(time).getFloats(0)
> 
> #Reflectivity of ground(sg)
> CSA = valuesCSA[0][crd]
> print "CSA: ", CSA, " knotlength: ", knotLength * 10
> if (CSA > (knotLength-1*10)):
> CSA = (knotLength-1)*10
> if CSA < 0:
> CSA = 0
> sg = function.value(Double(CSA))
> 
> singleB = currTimeSequence.getSample(time)
> lineValues  = singleB.getFloats(0)
> 
> #Reflectivity of current pixel
> sc = lineValues[0][crd]/100
> print "Reflectivity of current pixel: ", sc
> 
> # Cloud index(n)  rho_cloud =0.81 (Dagestad, personal comm.) or use 
> pth-percentile through user_pCloud
> #n = (lineValuesA[0][i] * 0.01  - 0.17) / (0.81-0.17) # 0.17 is from weather 
> station data
> n = (sc - sg) / (rc - sg)
> if n <= -0.2:
> print "n <= -0.2"
> K = 1.2
> print "n: ",n," K: ", K
> elif n > 1.1:
> print "n >= 1.1"
> K = 0.05
> print "n: ",n," K: ", K
> elif n > 0.80:
> print "n > 0.80"
> K = 2.0667 - 3.6667 * n + 1.6667 * n * n
> print "n: ",n," K: ", K
> else:
> print "n <= 0.80"
> K = 1 - n
> print "n: ",n," K: ", K
> 
> #Correction of k#new Kch = Kch ? 0.01(8 TST ? 104)
> #Kc = K - 0.01 * (8 * utc - 104)
> #if Kc > 1.2:
> #  Kcr = 1.2
> #elif Kc < 0.05:
> #  Kcr = 0.05
> #else:
> #  Kcr = Kc
> 
> try:
> lineValues[0][crd] = K
> except:
> lineValues[0][crd] = Float.NaN #error, fill pixel value with "missing" or NaN
> singleB.setSamples(lineValues)
> newd.setSample(time,singleB)
> return newd
> 
> 
> 


Ticket Details
===================
Ticket ID: VVW-585732
Department: Support IDV
Priority: Normal
Status: Open