Hi Tom,
The reason I implemented the "if temp==NaN...." logic is because IDV was
giving me 'NaN' errors. I can't remember exactly the reason, but I
inserted this code as a 'quick fix' - perhaps not good programming practice
but I was eager to see if anything was working so I could get some results
to work with and it has remained.
Below is my code for calculating the wet-bulb temperature at each point.
I'm happy for others to use it and perhaps modify it to come up with
something better.
It works by adiabatically and isobarically cooling the dry bulb temperature
(or starting from a first guess temperature) to the wetbulb temperature,
which is implicitly defined by its saturation mixing ratio. Perhaps there
is a better and faster way of doing this with IDV?
#Copyright (C) Paul Graham, December 2013. This library is free software;
#you can redistribute it and/or modify it under the terms of the GNU Lesser
#General Public License as published by the Free Software Foundation
def bolton(t):
"""Estimates the saturation vapour pressure in hPa at a given temperature,
T, in celcius"""
from java.lang import Math
# return saturation vapour pressure at a given temperature in celcius
es = 6.112*Math.exp(17.67*t/(t+243.5))
return es
def findwetbulb(t,td,p):
"""Iteratively finds the wet bulb temperature in celcius at pressure
level, p in hPa, given the temperature, t in celcius, and dew point, td in
celcius"""
R = 287
Cp = 1004
L = 2.5E6
r = 0.622*bolton(td)/p
tw = t
limit = r+(Cp/L)*t
guess = (Cp/L)*tw+(0.622/p)*bolton(tw)
if(str(float(limit))!='nan' and str(float(guess))!='nan'):
while(guess>limit):
tw=tw-0.1
guess = (Cp/L)*tw+(0.622/p)*bolton(tw)
else: tw = java.lang.Float.NaN
return tw
def kelvinToCelcius(tk):
"""Returns the temperature in celcius given a temperature, tk, in
kelvin"""
return tk - 273
Paul
On 11 December 2013 01:50, Tom Whittaker <whittaker@xxxxxxxx> wrote:
> Hi Paul --
>
> I'm not clear on the logic of "(if temp == NaN) then compute
> WetBulb"...or am I mis-reading this?
>
> Can you send me (off-line?) your web bulb computation formula? I'd
> like to see what additional logic is in there that prevents you from
> using VisAD "built-in" arithmetic. (I'm always looking for ideas for
> new functions to put into JPythonMethods).
>
> tom
>
> On Mon, Dec 9, 2013 at 3:59 PM, Paul Graham <meteorpaul@xxxxxxxxx> wrote:
> > Hi Tom,
> >
> > Thanks for your help. I figured out a method after some experimentation
> and
> > realised I could not use the visad logic. My algorithm works point-wise,
> > which requires having to step through each point individually.
> Essentially,
> > I modified existing code for substitution, found in IDV's Jython library:
> >
> > def wetbulb(temp,dewpt):
> > from java.lang import Math
> > tempData = temp.clone();
> > dewptData = dewpt.clone();
> > presData = extractPressureFromNWPGrid(tempData);
> > wetbulbData = temp.clone();
> > if (GridUtil.isTimeSequence(tempData)):
> > for t in range(tempData.getDomainSet().getLength()):
> > tempValues = tempData[t]
> > dewptValues = dewptData[t]
> > presValues = presData[t]
> > wetbulbRangeObj = wetbulbData.getSample(t)
> > newValues = wetbulbRangeObj.getFloats(0)
> > print "Processing wetbulb temperature for time ", t
> > for i in range(len(tempValues)):
> > if(str(float(tempValues[i].getValue()))!='nan'):
> > newValues[0][i] =
> >
> findwetbulb(kelvinToCelcius(tempValues[i].getValue()),kelvinToCelcius(dewptValues[i].getValue()),presData[i].getValue())
> > wetbulbRangeObj.setSamples(newValues,1)
> > wetbulbData = newUnit(noUnit(wetbulbData),'Wetbulb temperature in
> > celcius','C')
> > return wetbulbData
> >
> >
> > Paul
> >
> >
> >
> >
> > On 10 December 2013 01:12, Tom Whittaker <whittaker@xxxxxxxx> wrote:
> >>
> >> Hi Paul --
> >>
> >> If your function has no logic that needs to be applied to each point,
> >> then if the data were read in via the IDV they would be "VisAD Data
> >> objects" and you should just be able to do something like:
> >>
> >> a = b + c - d
> >>
> >> where "b", "c" and "d" are 3D Data objects (or maybe 4D if they
> >> include "time") and it will do the computation on every point. and put
> >> the result in "a" (creating a new Data object with characteristics
> >> like that of "b" (the first object in the computation).
> >>
> >> If you have logic, then you might have to iterate over every point...I
> >> say "might" because sometimes some of the built-in methods like "mask"
> >> can be used to simulate logic.
> >>
> >> Hope that helps.
> >>
> >> tom
> >>
> >>
> >>
> >> On Sun, Dec 8, 2013 at 9:39 PM, Paul Graham <meteorpaul@xxxxxxxxx>
> wrote:
> >> > Hi IDV'ers,
> >> >
> >> > I have written a Jython function which takes temperature, dewpoint and
> >> > pressure as arguments and evaluated at each point in a 3d domain to
> >> > calculate the wet bulb temperature. What is the best way using Jython
> >> > to
> >> > iterate over my domain so I can set the wetbulb result for each point?
> >> > Eg. wetbulbresult[i][j][k] =
> >> > findwetbulb(temperature[i][j][k],dewpoint[i][j][k],pressure[i][j][k]),
> >> > where i,j are subscripts for the latitude and longitude and k for the
> >> > height.
> >> >
> >> > Thanks in advance,
> >> >
> >> > Paul
> >> > _______________________________________________
> >> > idvusers mailing list
> >> > idvusers@xxxxxxxxxxxxxxxx
> >> > For list information, to unsubscribe, visit:
> >> > http://www.unidata.ucar.edu/mailing_lists/
> >>
> >>
> >>
> >> --
> >> Tom Whittaker
> >> University of Wisconsin-Madison
> >> Space Science & Engineering Center (SSEC)
> >> Cooperative Institute for Meteorological Satellite Studies (CIMSS)
> >> 1225 W. Dayton Street
> >> Madison, WI 53706 USA
> >> ph: +1 608 262 2759
> >
> >
>
>
>
> --
> Tom Whittaker
> University of Wisconsin-Madison
> Space Science & Engineering Center (SSEC)
> Cooperative Institute for Meteorological Satellite Studies (CIMSS)
> 1225 W. Dayton Street
> Madison, WI 53706 USA
> ph: +1 608 262 2759
>
--
Paul Graham
Email: meteorpaul@xxxxxxxxx
Mobile: 0403003784
