uWxUtils
uWxUtils.zip contains Delphi Object Pascal source code files containing
all the weather algorithms used in VpLive and VpPressureCalc. The most
important algorithms convert pressure between sensor pressure, station
pressure, altimeter, and sea level reduced pressure. There are three
files: uWxUtils.pas is the main file containing the various functions,
all of which take and return values in metric units, uWxUtilsUS.pas
wraps the functions in uWxUtils.pas with functions that take and return
values in US units, and uWxUtilsVP.pas contains pressure conversion
functions specific to Vantage Pro weather stations. The
uWxUtils.pas file contents are included at the bottom
of this page, with hyperlinked references. This code may be used
by anyone finding it useful, including commercial use.
uWxUtils most recent version
June 2006 - This code was previously called uWeatherUtils, and consisted of a single
file. You can download the previous version
here. The newer version was
reorganized to make a clearer division between generic weather functions,
and those specific to the Vantage Pro. I also wanted to make a clearer
division between metric and US unit functions, and to put all the core
functions in the metric library. I've also made some minor changes,
and added functions for computing Heat Index and Humidex.
Nov. 2007 - Broken hyperlinks to WMO website were fixed. An error in the
wind chill formula (thanks to Juergen for pointing this out) was corrected.
The formula was using the value 13.112, but the correct value is 13.12.
The magnitude of the error is less than a degree. I also added links to
sources for this forumula. Also corrected a bug in the HumidityCorrection
function. It was always using the vaDavis vapor pressure algorithm instead
of using the algorithm passed in.
(pas to html formatting via Peter Johnson's
tool at DelphiDabbler
unit uWxUtils;
interface
type
TWxReal = double;
TWxHumidity = byte;
TSLPAlgorithm = (
paDavisVP,
paUnivie,
paManBar
);
TAltimeterAlgorithm = (
aaASOS,
aaASOS2,
aaMADIS,
aaNOAA,
aaWOB,
aaSMT
);
TVapAlgorithm = (
vaDavisVp,
vaBuck,
vaBuck81,
vaBolton,
vaTetenNWS,
vaTetenMurray,
vaTeten
);
const
DefaultSLPAlgorithm = paManBar;
DefaultAltimeterAlgorithm = aaMADIS;
DefaultVapAlgorithm = vaBolton;
type
TWxUtils = class(TObject)
private
public
class function StationToSensorPressure(pressureHPa: TWxReal;
sensorElevationM: TWxReal; stationElevationM: TWxReal;
currentTempC: TWxReal): TWxReal;
class function StationToAltimeter(pressureHPa: TWxReal; elevationM: TWxReal;
algorithm: TAltimeterAlgorithm = DefaultAltimeterAlgorithm): TWxReal;
class function StationToSeaLevelPressure(pressureHPa: TWxReal; elevationM: TWxReal;
currentTempC: TWxReal; meanTempC: TWxReal; humidity: TWxHumidity;
algorithm: TSLPAlgorithm = DefaultSLPAlgorithm): TWxReal;
class function SensorToStationPressure(pressureHPa: TWxReal;
sensorElevationM: TWxReal; stationElevationM: TWxReal;
currentTempC: TWxReal): TWxReal;
class function SeaLevelToStationPressure(pressureHPa: TWxReal; elevationM: TWxReal;
currentTempC: TWxReal; meanTempC: TWxReal; humidity: TWxHumidity;
algorithm: TSLPAlgorithm = DefaultSLPAlgorithm): TWxReal;
class function PressureReductionRatio(pressureHPa: TWxReal; elevationM: TWxReal;
currentTempC: TWxReal; meanTempC: TWxReal; humidity: TWxHumidity;
algorithm: TSLPAlgorithm = DefaultSLPAlgorithm): TWxReal;
class function ActualVaporPressure(tempC: TWxReal; humidity: TWxHumidity;
algorithm: TVapAlgorithm = DefaultVapAlgorithm): TWxReal;
class function SaturationVaporPressure(tempC: TWxReal;
algorithm: TVapAlgorithm = DefaultVapAlgorithm): TWxReal;
class function MixingRatio(pressureHPa: TWxReal; tempC: TWxReal; humidity: TWxHumidity): TWxReal;
class function VirtualTempK(pressureHPa: TWxReal; tempC: TWxReal; humidity: TWxHumidity): TWxReal;
class function HumidityCorrection(tempC: TWxReal; elevationM: TWxReal; humidity: TWxHumidity;
algorithm: TVapAlgorithm = DefaultVapAlgorithm): TWxReal;
class function DewPoint(tempC: TWxReal; humidity: TWxHumidity;
algorithm: TVapAlgorithm = DefaultVapAlgorithm): TWxReal;
class function WindChill(tempC: TWxReal; windSpeedKmph: TWxReal): TWxReal;
class function HeatIndex(tempC: TWxReal; humidity: TWxHumidity): TWxReal;
class function Humidex(tempC: TWxReal; humidity: TWxHumidity): TWxReal;
class function GeopotentialAltitude(geometricAltitudeM: TWxReal): TWxReal;
class function FToC(value: TWxReal): TWxReal;
class function CToF(value: TWxReal): TWxReal;
class function CToK(value: TWxReal): TWxReal;
class function KToC(value: TWxReal): TWxReal;
class function FToR(value: TWxReal): TWxReal;
class function RToF(value: TWxReal): TWxReal;
class function InToHPa(value: TWxReal): TWxReal;
class function HPaToIn(value: TWxReal): TWxReal;
class function FtToM(value: TWxReal): TWxReal;
class function MToFt(value: TWxReal): TWxReal;
class function InToMm(value: TWxReal): TWxReal;
class function MmToIn(value: TWxReal): TWxReal;
class function MToKm(value: TWxReal): TWxReal;
class function KmToM(value: TWxReal): TWxReal;
class function Power(const Base, Exponent: TWxReal): TWxReal;
class function Power10(const exponent: TWxReal): TWxReal;
end;
implementation
const
gravity = 9.80665;
uGC = 8.31432;
moleAir = 0.0289644;
moleWater = 0.01801528;
gasConstantAir = uGC/moleAir;
standardSLP = 1013.25;
standardSlpInHg = 29.921;
standardTempK = 288.15;
earthRadius45 = 6356.766;
standardLapseRate = 0.0065;
standardLapseRateFt = standardLapseRate * 0.3048;
vpLapseRateUS = 0.00275;
manBarLapseRate = 0.0117;
class function TWxUtils.StationToSensorPressure(pressureHPa: TWxReal;
sensorElevationM: TWxReal; stationElevationM: TWxReal;
currentTempC: TWxReal): TWxReal;
begin
Result := InToHPa(HPaToIn(pressureHPa) / Power10(0.00813 * MToFt(sensorElevationM - stationElevationM)
/ FToR(CToF(currentTempC))));
end;
class function TWxUtils.StationToAltimeter(PressureHPa: TWxReal; elevationM: TWxReal;
algorithm: TAltimeterAlgorithm = DefaultAltimeterAlgorithm): TWxReal;
var
geopEl: TWxReal;
k1, k2: TWxReal;
begin
case algorithm of
aaASOS:
begin
Result := InToHPa(Power(Power(HPaToIn(pressureHPa), 0.1903) + (1.313E-5 * MToFt(elevationM)), 5.255));
end;
aaASOS2:
begin
geopEl := GeopotentialAltitude(elevationM);
k1 := standardLapseRate * gasConstantAir / gravity;
k2 := 8.41728638E-5;
Result := Power(Power(pressureHPa, k1) + (k2 * geopEl), 1/k1);
end;
aaMADIS:
begin
k1 := 0.190284;
k2 := 8.4184960528E-5;
Result := Power(Power(pressureHPa - 0.3, k1) + (k2 * elevationM), 1/k1);
end;
aaNOAA:
begin
k1 := 0.190284;
k2 := 8.42288069E-5;
Result := (pressureHPa - 0.3) * Power(1 + (k2 * (elevationM / Power(pressureHPa - 0.3, k1))), 1/k1);
end;
aaWOB:
begin
k1 := standardLapseRate * gasConstantAir / gravity;
k2 := 1.312603E-5;
Result := InToHPa(Power(Power(HPaToIn(pressureHPa), k1) + (k2 * MToFt(elevationM)), 1/k1));
end;
aaSMT:
begin
k1 := 0.190284;
k2 := 4.30899E-5;
geopEl := GeopotentialAltitude(elevationM);
Result := InToHPa((HPaToIn(pressureHPa) - 0.01)
* Power(1 + (k2 * (geopEl / Power(HPaToIn(pressureHPa) - 0.01, k1))), 1/k1));
end;
else Result := pressureHPa;
end;
end;
class function TWxUtils.StationToSeaLevelPressure(pressureHPa: TWxReal; elevationM: TWxReal;
currentTempC: TWxReal; meanTempC: TWxReal; humidity: TWxHumidity;
algorithm: TSLPAlgorithm = DefaultSLPAlgorithm): TWxReal;
begin
Result := pressureHPa * PressureReductionRatio(pressureHPa, elevationM, currentTempC,
meanTempC, humidity, algorithm);
end;
class function TWxUtils.SensorToStationPressure(pressureHPa: TWxReal;
sensorElevationM: TWxReal; stationElevationM: TWxReal;
currentTempC: TWxReal): TWxReal;
begin
Result := InToHPa(HPaToIn(pressureHPa) * Power10(0.00813 * MToFt(sensorElevationM - stationElevationM)
/ FToR(CToF(currentTempC))));
end;
class function TWxUtils.SeaLevelToStationPressure(pressureHPa: TWxReal; elevationM: TWxReal;
currentTempC: TWxReal; meanTempC: TWxReal; humidity: TWxHumidity;
algorithm: TSLPAlgorithm = DefaultSLPAlgorithm): TWxReal;
begin
Result := pressureHPa / PressureReductionRatio(pressureHPa, elevationM, currentTempC,
meanTempC, humidity, algorithm);
end;
class function TWxUtils.PressureReductionRatio(pressureHPa: TWxReal; elevationM: TWxReal;
currentTempC: TWxReal; meanTempC: TWxReal; humidity: TWxHumidity;
algorithm: TSLPAlgorithm = DefaultSLPAlgorithm): TWxReal;
var
geopElevationM: TWxReal;
hCorr: TWxReal;
begin
case algorithm of
paUnivie:
begin
geopElevationM := GeopotentialAltitude(elevationM);
Result := Exp(((gravity/gasConstantAir) * geopElevationM)
/ (VirtualTempK(pressureHPa, meanTempC, humidity) + (geopElevationM * standardLapseRate/2)));
end;
paDavisVp:
begin
if (humidity > 0) then begin
hCorr := (9/5) * HumidityCorrection(currentTempC, elevationM, humidity, vaDavisVP);
end else begin
hCorr := 0;
end;
Result := Power(10, (MToFt(elevationM) / (122.8943111 * (CToF(meanTempC) + 460 + (MToFt(elevationM) * vpLapseRateUS/2) + hCorr))));
end;
paManBar:
begin
if (humidity > 0) then begin
hCorr := (9/5) * HumidityCorrection(currentTempC, elevationM, humidity, vaBuck);
end else begin
hCorr := 0;
end;
geopElevationM := GeopotentialAltitude(elevationM);
Result := Exp(geopElevationM * 6.1454E-2 / (CToF(meanTempC) + 459.7 + (geopElevationM * manBarLapseRate / 2) + hCorr));
end;
else Result := 1;
end;
end;
class function TWxUtils.ActualVaporPressure(tempC: TWxReal; humidity: TWxHumidity;
algorithm: TVapAlgorithm = DefaultVapAlgorithm): TWxReal;
begin
result := (humidity * SaturationVaporPressure(tempC, algorithm)) / 100;
end;
class function TWxUtils.SaturationVaporPressure(tempC: TWxReal;
algorithm: TVapAlgorithm = DefaultVapAlgorithm): TWxReal;
begin
case algorithm of
vaDavisVp: Result := 6.112 * Exp((17.62 * tempC)/(243.12 + tempC));
vaBuck: Result := 6.1121 * Exp((18.678 - (tempC/234.5)) * tempC / (257.14 + tempC));
vaBuck81: Result := 6.1121 * Exp((17.502 * tempC)/(240.97 + tempC));
vaBolton: Result := 6.112 * Exp(17.67 * tempC / (tempC + 243.5));
vaTetenNWS: Result := 6.112 * Power(10,(7.5 * tempC / (tempC + 237.7)));
vaTetenMurray: Result := Power(10, (7.5 * tempC / (237.5 + tempC)) + 0.7858);
vaTeten: Result := 6.1078 * Power(10, (7.5 * tempC / (tempC + 237.3)));
else Result := 0;
end;
end;
class function TWxUtils.MixingRatio(pressureHPa: TWxReal; tempC: TWxReal;
humidity: TWxHumidity): TWxReal;
var
vapPres: TWxReal;
const
k1 = moleWater / moleAir;
begin
vapPres := ActualVaporPressure(tempC, humidity, vaBuck);
Result := 1000 * ((k1 * vapPres) / (pressureHPa - vapPres));
end;
class function TWxUtils.VirtualTempK(pressureHPa: TWxReal; tempC: TWxReal;
humidity: TWxHumidity): TWxReal;
var
vapPres: TWxReal;
const
epsilon = 1 - (moleWater / moleAir);
begin
vapPres := ActualVaporPressure(tempC, humidity, vaBuck);
Result := (CtoK(tempC)) / (1-(epsilon * (vapPres/pressureHPa)));
end;
class function TWxUtils.HumidityCorrection(tempC: TWxReal; elevationM: TWxReal; humidity: TWxHumidity;
algorithm: TVapAlgorithm = DefaultVapAlgorithm): TWxReal;
var
vapPress: TWxReal;
begin
vapPress := ActualVaporPressure(tempC, humidity, algorithm);
Result := (vappress * ((2.8322E-9 * Sqr(elevationM)) + (2.225E-5 * elevationM) + 0.10743));
end;
class function TWxUtils.DewPoint(tempC: TWxReal; humidity: TWxHumidity;
algorithm: TVapAlgorithm = DefaultVapAlgorithm): TWxReal;
var
lnVapor: TWxReal;
begin
lnVapor := Ln(ActualVaporPressure(tempc, humidity, algorithm));
case algorithm of
vaDavisVp: Result := ((243.12 * LnVapor) - 440.1) / (19.43 - LnVapor);
else
Result := ((237.7 * LnVapor) - 430.22) / (19.08 - LnVapor);
end;
end;
class function TWxUtils.WindChill(tempC: TWxReal; windSpeedKmph: TWxReal): TWxReal;
var
windPow: TWxReal;
begin
if ((tempC >= 10.0) or (windSpeedKmph <= 4.8)) then begin
Result := tempC;
end else begin
windPow := Power(windSpeedKmph, 0.16);
Result := 13.12 + (0.6215 * tempC) - (11.37 * windPow) + (0.3965 * tempC * windPow);
end;
if (Result > tempC) then Result := tempC;
end;
class function TWxUtils.HeatIndex(tempC: TWxReal; humidity: TWxHumidity): TWxReal;
var
tempF: TWxReal;
tSqrd: TWxReal;
hum: TWxReal;
hSqrd: TWxReal;
begin
tempF := CToF(tempC);
if (tempF < 80) then begin
Result := tempF;
end else begin
tSqrd := tempF * tempF;
hum := humidity;
hSqrd := hum * hum;
Result := 0 - 42.379 + (2.04901523 * tempF) + (10.14333127 * humidity)
- (0.22475541 * tempF * humidity) - (0.00683783 * tSqrd)
- (0.05481717 * hSqrd) + (0.00122874 * tSqrd * humidity)
+ (0.00085282 * tempF * hSqrd) - (0.00000199 * tSqrd * hSqrd);
if ((humidity < 13) and (tempF >= 80) and (tempF <= 112)) then begin
Result := Result - ((13 - humidity)/4) * Sqrt((17 - Abs(tempf - 95))/17);
end else if ((humidity > 85) and (tempF >= 80) and (tempF <= 87)) then begin
Result := Result + ((humidity - 85)/10) * ((87 - tempF)/5);
end;
end;
Result := FToC(Result);
end;
class function TWxUtils.Humidex(tempC: TWxReal; humidity: TWxHumidity): TWxReal;
begin
Result := tempC + ((5/9) * (ActualVaporPressure(tempC, humidity, vaTetenNWS) - 10.0));
end;
class function TWxUtils.GeopotentialAltitude(geometricAltitudeM: TWxReal): TWxReal;
begin
Result := (earthRadius45 * 1000 * geometricAltitudeM) / ((earthRadius45 * 1000) + geometricAltitudeM);
end;
class function TWxUtils.FtoC(value: TWxReal): TWxReal;
begin
Result := ((value - 32) * 5) / 9;
end;
class function TWxUtils.CtoF(value: TWxReal): TWxReal;
begin
Result := ((value * 9) / 5) + 32;
end;
class function TWxUtils.CtoK(value: TWxReal): TWxReal;
begin
Result := 273.15 + value;
end;
class function TWxUtils.KtoC(value: TWxReal): TWxReal;
begin
Result := value - 273.15;
end;
class function TWxUtils.FToR(value: TWxReal): TWxReal;
begin
Result := value + 459.67;
end;
class function TWxUtils.RToF(value: TWxReal): TWxReal;
begin
Result := value - 459.67;
end;
class function TWxUtils.InToHPa(value: TWxReal): TWxReal;
begin
Result := value / 0.02953;
end;
class function TWxUtils.HPaToIn(value: TWxReal): TWxReal;
begin
Result := value * 0.02953;
end;
class function TWxUtils.FtToM(value: TWxReal): TWxReal;
begin
Result := value * 0.3048;
end;
class function TWxUtils.MToFt(value: TWxReal): TWxReal;
begin
Result := value / 0.3048;
end;
class function TWxUtils.InToMm(value: TWxReal): TWxReal;
begin
Result := value * 25.4;
end;
class function TWxUtils.MmToIn(value: TWxReal): TWxReal;
begin
Result := value / 25.4;
end;
class function TWxUtils.MToKm(value: TWxReal): TWxReal;
begin
Result := value * 1.609344;
end;
class function TWxUtils.KmToM(value: TWxReal): TWxReal;
begin
Result := value / 1.609344;
end;
class function TWxUtils.Power(const base, exponent: TWxReal): TWxReal;
begin
if exponent = 0.0 then begin
Result := 1.0;
end else if (base = 0.0) and (exponent > 0.0) then begin
Result := 0.0;
end else begin
Result := Exp(exponent * Ln(base));
end;
end;
class function TWxUtils.Power10(const exponent: TWxReal): TWxReal;
const
ln10 = 2.302585093;
begin
if (exponent = 0.0) then begin
Result := 1.0;
end else begin
Result := Exp(exponent * ln10);
end;
end;
end.