Pressure to altitude conversion fixed (#2483)

Author: pgrawehr

src/Iot.Device.Bindings/CompatibilitySuppressions.xml

@@ -1,6 +1,13 @@
 <?xml version="1.0" encoding="utf-8"?>
 <!-- https://learn.microsoft.com/dotnet/fundamentals/package-validation/diagnostic-ids -->
 <Suppressions xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:xsd="http://www.w3.org/2001/XMLSchema">
+  <Suppression>
+    <DiagnosticId>CP0002</DiagnosticId>
+    <Target>F:Iot.Device.Common.WeatherHelper.MeanSeaLevel</Target>
+    <Left>lib/net8.0/Iot.Device.Bindings.dll</Left>
+    <Right>lib/net8.0/Iot.Device.Bindings.dll</Right>
+    <IsBaselineSuppression>true</IsBaselineSuppression>
+  </Suppression>
   <Suppression>
     <DiagnosticId>CP0002</DiagnosticId>
     <Target>F:Iot.Device.Tca955x.Tca955x.DefaultI2cAdress</Target>
@@ -15,6 +22,13 @@
     <Right>lib/net8.0/Iot.Device.Bindings.dll</Right>
     <IsBaselineSuppression>true</IsBaselineSuppression>
   </Suppression>
+  <Suppression>
+    <DiagnosticId>CP0002</DiagnosticId>
+    <Target>M:Iot.Device.Common.WeatherHelper.CalculateAltitude(UnitsNet.Pressure,UnitsNet.Temperature)</Target>
+    <Left>lib/net8.0/Iot.Device.Bindings.dll</Left>
+    <Right>lib/net8.0/Iot.Device.Bindings.dll</Right>
+    <IsBaselineSuppression>true</IsBaselineSuppression>
+  </Suppression>
   <Suppression>
     <DiagnosticId>CP0002</DiagnosticId>
     <Target>M:Iot.Device.Nmea0183.Ais.TrackEstimationParameters.get_WarningDistance</Target>

src/devices/Bmp180/Bmp180.cs

@@ -98,7 +98,7 @@ public Pressure ReadPressure()
         /// <returns>
         ///  Height in meters above sea level
         /// </returns>
-        public Length ReadAltitude() => ReadAltitude(WeatherHelper.MeanSeaLevel);
+        public Length ReadAltitude() => ReadAltitude(WeatherHelper.MeanSeaLevelPressure);
 
         /// <summary>
         ///  Calculates the pressure at sea level when given a known altitude

src/devices/Bmp180/samples/Program.cs

@@ -27,8 +27,8 @@
 Console.WriteLine($"Pressure: {preValue.Hectopascals:0.##}hPa");
 
 // Note that if you already have the pressure value and the temperature, you could also calculate altitude by
-// calling WeatherHelper.CalculateAltitude(preValue, Pressure.MeanSeaLevel, tempValue) which would be more performant.
-Length altValue = i2cBmp280.ReadAltitude(WeatherHelper.MeanSeaLevel);
+// calling WeatherHelper.CalculateAltitude(preValue, WeatherHelper.MeanSeaLevelPressure, tempValue) which would be more performant.
+Length altValue = i2cBmp280.ReadAltitude(WeatherHelper.MeanSeaLevelPressure);
 
 Console.WriteLine($"Altitude: {altValue:0.##}m");
 Thread.Sleep(1000);
@@ -43,6 +43,6 @@
 Console.WriteLine($"Pressure: {preValue.Hectopascals:0.##}hPa");
 
 // Note that if you already have the pressure value and the temperature, you could also calculate altitude by
-// calling WeatherHelper.CalculateAltitude(preValue, Pressure.MeanSeaLevel, tempValue) which would be more performant.
-altValue = i2cBmp280.ReadAltitude(WeatherHelper.MeanSeaLevel);
+// calling WeatherHelper.CalculateAltitude(preValue, WeatherHelper.MeanSeaLevelPressure, tempValue) which would be more performant.
+altValue = i2cBmp280.ReadAltitude(WeatherHelper.MeanSeaLevelPressure);
 Console.WriteLine($"Altitude: {altValue:0.##}m");

src/devices/Bmxx80/Bmx280Base.cs

@@ -176,7 +176,7 @@ public bool TryReadAltitude(Pressure seaLevelPressure, out Length altitude)
         /// Contains <see cref="double.NaN"/> otherwise.
         /// </param>
         /// <returns><code>true</code> if pressure measurement was not skipped, otherwise <code>false</code>.</returns>
-        public bool TryReadAltitude(out Length altitude) => TryReadAltitude(WeatherHelper.MeanSeaLevel, out altitude);
+        public bool TryReadAltitude(out Length altitude) => TryReadAltitude(WeatherHelper.MeanSeaLevelPressure, out altitude);
 
         /// <summary>
         /// Get the current status of the device.

src/devices/Bmxx80/samples/Bme280/Bme280.sample.cs

@@ -15,7 +15,7 @@
 // bus id on the raspberry pi 3
 const int busId = 1;
 // set this to the current sea level pressure in the area for correct altitude readings
-Pressure defaultSeaLevelPressure = WeatherHelper.MeanSeaLevel;
+Pressure defaultSeaLevelPressure = WeatherHelper.MeanSeaLevelPressure;
 
 I2cConnectionSettings i2cSettings = new(busId, Bme280.DefaultI2cAddress);
 using I2cDevice i2cDevice = I2cDevice.Create(i2cSettings);

src/devices/Bmxx80/samples/Bme680/Bme680.sample.cs

@@ -13,7 +13,7 @@
 // The I2C bus ID on the Raspberry Pi 3.
 const int busId = 1;
 // set this to the current sea level pressure in the area for correct altitude readings
-Pressure defaultSeaLevelPressure = WeatherHelper.MeanSeaLevel;
+Pressure defaultSeaLevelPressure = WeatherHelper.MeanSeaLevelPressure;
 
 I2cConnectionSettings i2cSettings = new(busId, Bme680.DefaultI2cAddress);
 I2cDevice i2cDevice = I2cDevice.Create(i2cSettings);

src/devices/Bmxx80/samples/Bmp280/Bmp280.sample.cs

@@ -16,7 +16,7 @@
 // bus id on the raspberry pi 3 and 4
 const int busId = 1;
 // set this to the current sea level pressure in the area for correct altitude readings
-Pressure defaultSeaLevelPressure = WeatherHelper.MeanSeaLevel;
+Pressure defaultSeaLevelPressure = WeatherHelper.MeanSeaLevelPressure;
 
 I2cConnectionSettings i2cSettings = new(busId, Bmp280.DefaultI2cAddress);
 I2cDevice i2cDevice = I2cDevice.Create(i2cSettings);

src/devices/Common/Iot/Device/Common/WeatherHelper.cs

@@ -15,22 +15,30 @@ public static class WeatherHelper
         /// <summary>
         /// Gas constant of dry Air, J / (kg * K)
         /// </summary>
-        internal const double SpecificGasConstantOfAir = 287.058;
+        public const double SpecificGasConstantOfAir = 287.058;
 
         /// <summary>
         /// Gas constant of vapor, J / (kg * K)
         /// </summary>
-        internal const double SpecificGasConstantOfVapor = 461.523;
+        public const double SpecificGasConstantOfVapor = 461.523;
 
         /// <summary>
         /// Default atmospheric temperature gradient = 0.0065K/m (or 0.65K per 100m)
         /// </summary>
-        internal const double DefaultTemperatureGradient = 0.0065;
+        public const double DefaultTemperatureGradient = 0.0065;
+
+        /// <summary>
+        /// The default assumed temperature at sea level.
+        /// This is used for <see cref="CalculateAltitude(Pressure, Pressure)"/> since it is not possible to calculate the temperature at sea level
+        /// for this formula, as we do not know the altitude (it is the result of that computation, after all). The QNH value does consider the
+        /// temperature, though.
+        /// </summary>
+        public static Temperature DefaultSeaLevelTemperature => Temperature.FromDegreesCelsius(15);
 
         /// <summary>
         /// The mean sea-level pressure (MSLP) is the average atmospheric pressure at mean sea level
         /// </summary>
-        public static readonly Pressure MeanSeaLevel = Pressure.FromPascals(101325);
+        public static Pressure MeanSeaLevelPressure => Pressure.FromPascals(101325);
 
         #region TemperatureAndRelativeHumidity
 
@@ -191,47 +199,46 @@ public static RelativeHumidity GetRelativeHumidityFromActualAirTemperature(Tempe
         #endregion TemperatureAndRelativeHumidity
 
         #region Pressure
-        // Formula  from https://de.wikipedia.org/wiki/Barometrische_Höhenformel#Internationale_Höhenformel, solved
-        // for different parameters
 
         /// <summary>
-        /// Calculates the altitude in meters from the given pressure, sea-level pressure and air temperature.
+        /// Calculates the altitude in meters from the given pressure and sea-level pressure.
         /// </summary>
         /// <param name="pressure">The pressure at the point for which altitude is being calculated</param>
-        /// <param name="seaLevelPressure">The sea-level pressure</param>
-        /// <param name="airTemperature">The dry air temperature at the point for which altitude is being calculated</param>
-        /// <returns>The altitude</returns>
-        public static Length CalculateAltitude(Pressure pressure, Pressure seaLevelPressure, Temperature airTemperature)
+        /// <param name="seaLevelPressure">The sea-level pressure. In aviation, this is called the QNH value and provided by weather reports or ATC</param>
+        /// <remarks>
+        /// Formula  from https://de.wikipedia.org/wiki/Barometrische_Höhenformel#Internationale_Höhenformel, solved
+        /// </remarks>
+        /// <returns>The altitude over mean sea level</returns>
+        public static Length CalculateAltitude(Pressure pressure, Pressure seaLevelPressure)
         {
-            double meters = ((Math.Pow(seaLevelPressure.Pascals / pressure.Pascals, 1 / 5.255) - 1) * airTemperature.Kelvins) / DefaultTemperatureGradient;
+            double meters = (Math.Pow(seaLevelPressure.Pascals / pressure.Pascals, -1.0 / 5.255) - 1) * (DefaultSeaLevelTemperature.Kelvins / -DefaultTemperatureGradient);
             return Length.FromMeters(meters);
         }
 
         /// <summary>
-        /// Calculates the altitude in meters from the given pressure and air temperature. Assumes mean sea-level pressure.
-        /// </summary>
-        /// <param name="pressure">The pressure at the point for which altitude is being calculated</param>
-        /// <param name="airTemperature">The dry air temperature at the point for which altitude is being calculated</param>
-        /// <returns>The altitude</returns>
-        public static Length CalculateAltitude(Pressure pressure, Temperature airTemperature) =>
-             CalculateAltitude(pressure, MeanSeaLevel, airTemperature);
-
-        /// <summary>
-        /// Calculates the altitude in meters from the given pressure and sea-level pressure. Assumes temperature of 15C.
+        /// Calculates the altitude in meters from the given pressure, sea-level pressure and air temperature.
         /// </summary>
         /// <param name="pressure">The pressure at the point for which altitude is being calculated</param>
-        /// <param name="seaLevelPressure">The sea-level pressure</param>
-        /// <returns>The altitude</returns>
-        public static Length CalculateAltitude(Pressure pressure, Pressure seaLevelPressure) =>
-            CalculateAltitude(pressure, seaLevelPressure, Temperature.FromDegreesCelsius(15));
+        /// <param name="seaLevelPressure">The sea-level pressure. In aviation, this is called the QNH value and provided by weather reports or ATC</param>
+        /// <param name="temperatureAtObservation">Temperature at observation point. Since this depends on the altitude (which we don't know), we use an iterative approach here</param>
+        /// <returns>The altitude over mean sea level</returns>
+        public static Length CalculateAltitude(Pressure pressure, Pressure seaLevelPressure, Temperature temperatureAtObservation)
+        {
+            // The first iteration is to get an approximation of the current height, so we can do the temperature compensation
+            double meters = (Math.Pow(seaLevelPressure.Pascals / pressure.Pascals, -1.0 / 5.255) - 1) * (DefaultSeaLevelTemperature.Kelvins / -DefaultTemperatureGradient);
+            // It gets hotter when going down (when in the Troposphere)
+            Temperature temperatureAtSeaLevel = temperatureAtObservation + TemperatureDelta.FromDegreesCelsius(DefaultTemperatureGradient * meters);
+            double meters2 = (Math.Pow(seaLevelPressure.Pascals / pressure.Pascals, -1.0 / 5.255) - 1) * (temperatureAtSeaLevel.Kelvins / -DefaultTemperatureGradient);
+            return Length.FromMeters(meters2);
+        }
 
         /// <summary>
-        /// Calculates the altitude in meters from the given pressure. Assumes mean sea-level pressure and temperature of 15C.
+        /// Calculates the altitude in meters from the given pressure. Assumes mean sea-level pressure.
         /// </summary>
         /// <param name="pressure">The pressure at the point for which altitude is being calculated</param>
         /// <returns>The altitude</returns>
         public static Length CalculateAltitude(Pressure pressure) =>
-            CalculateAltitude(pressure, MeanSeaLevel, Temperature.FromDegreesCelsius(15));
+            CalculateAltitude(pressure, MeanSeaLevelPressure);
 
         /// <summary>
         /// Calculates the approximate sea-level pressure from given absolute pressure, altitude and air temperature.

src/devices/Common/tests/WeatherTests.cs

@@ -81,33 +81,29 @@ public void AbsoluteHumidityIsCalculatedCorrectly(double expected, double fahren
         }
 
         [Theory]
-        [InlineData(1011.22, 900)]
-        [InlineData(111.18, 1000)]
-        [InlineData(547.1, 950)]
+        [InlineData(988.50, 900)]
+        [InlineData(110.88, 1000)]
         public void AltitudeIsCalculatedCorrectlyAtMslpAndDefaultTemp(double expected, double hpa)
         {
             Length altitude = WeatherHelper.CalculateAltitude(Pressure.FromHectopascals(hpa));
-            Assert.Equal(expected, Math.Round(altitude.Meters, 2));
+            Assert.Equal(altitude.Meters, expected, 0.5);
         }
 
         [Theory]
-        [InlineData(1011.22, 900, 1013.25)]
-        [InlineData(111.18, 1000, 1013.25)]
-        [InlineData(547.1, 950, 1013.25)]
-        public void AltitudeIsCalculatedCorrectlyAtDefaultTemp(double expected, double hpa, double seaLevelHpa)
+        [InlineData(0, 1013.25, 1013.25, 15)]
+        [InlineData(1000.13, 898.75, 1013.25, 15)]
+        [InlineData(3000.43, 701.09, 1013.25, 15)]
+        [InlineData(540.43, 950, 1013.25, 15)]
+        public void AltitudeIsCalculatedCorrectly(double expected, double hpa, double seaLevelHpa, double celsius)
         {
             Length altitude = WeatherHelper.CalculateAltitude(Pressure.FromHectopascals(hpa), Pressure.FromHectopascals(seaLevelHpa));
-            Assert.Equal(expected, Math.Round(altitude.Meters, 2));
-        }
 
-        [Theory]
-        [InlineData(1011.22, 900, 1013.25, 15)]
-        [InlineData(111.18, 1000, 1013.25, 15)]
-        [InlineData(547.1, 950, 1013.25, 15)]
-        public void AltitudeIsCalculatedCorrectly(double expected, double hpa, double seaLevelHpa, double celsius)
-        {
-            Length altitude = WeatherHelper.CalculateAltitude(Pressure.FromHectopascals(hpa), Pressure.FromHectopascals(seaLevelHpa), Temperature.FromDegreesCelsius(celsius));
+            // These two formulas should be equivalent.
+            double meters1 = (Math.Pow(seaLevelHpa / hpa, -1 / 5.255) - 1) * ((celsius + 273.15) / -WeatherHelper.DefaultTemperatureGradient);
+            double meters2 = ((celsius + 273.15) / WeatherHelper.DefaultTemperatureGradient) * (1 - Math.Pow(hpa / seaLevelHpa, 1.0 / 5.255));
             Assert.Equal(expected, Math.Round(altitude.Meters, 2));
+            Assert.Equal(meters1, expected, 0.01);
+            Assert.Equal(meters2, expected, 0.01);
         }
 
         [Theory]
@@ -159,6 +155,17 @@ public void CalculateBarometricPressure(double measuredValue, double temperature
             Assert.Equal(expected, result.Hectopascals, 2);
         }
 
+        [Theory]
+        [InlineData(1113.25, -801.147)]
+        [InlineData(1013.25, 0.0)]
+        [InlineData(913.25, 867.955)]
+        [InlineData(813.25, 1816.617)]
+        public void Pressure2Msl(double pressureHPa, double expectedHeightMsl)
+        {
+            var heightMsl = WeatherHelper.CalculateAltitude(Pressure.FromHectopascals(pressureHPa));
+            Assert.Equal(heightMsl.Meters, expectedHeightMsl, 0.1);
+        }
+
         [Theory]
         [InlineData(950, 15, 546.89, 10, 1013.19)]
         [InlineData(950, 15, 546.89, 50, 1013.01)]
@@ -235,5 +242,23 @@ public void CalculateWindForce(double temperature, double windSpeed, double pres
             var density = WeatherHelper.CalculateWindForce(airDensity, Speed.FromKilometersPerHour(windSpeed), 1.0);
             Assert.Equal(expected, density.NewtonsPerSquareMeter, 1);
         }
+
+        [Theory]
+        [InlineData(0, 1013.25, 15, 0)]
+        [InlineData(0, 1020, 15, 0)]
+        [InlineData(0, 1020, -20, 0)]
+        [InlineData(400, 970, 15, 397.989)]
+        [InlineData(400, 970, -5, 399.047)]
+        [InlineData(2000, 940, -5, 1995.953)]
+        public void RoundTripAltitudeAndPressure(double originalAltitude, double originalPressure, double originalTemperature, double expectedResult)
+        {
+            // The expectedResult and the originalAltitude should ideally be equal.
+            Pressure qnh = WeatherHelper.CalculateBarometricPressure(Pressure.FromHectopascals(originalPressure),
+                Temperature.FromDegreesCelsius(originalTemperature), Length.FromMeters(originalAltitude));
+
+            Length altitudeResult = WeatherHelper.CalculateAltitude(Pressure.FromHectopascals(originalPressure), qnh, Temperature.FromDegreesCelsius(originalTemperature));
+
+            Assert.Equal(expectedResult, altitudeResult.Meters, 1E-3);
+        }
     }
 }

src/devices/Ft232H/samples/Program.cs

@@ -115,7 +115,7 @@ void I2cScan(Ftx232HDevice ft232h)
 void TestI2c(Ftx232HDevice ft232h)
 {
     // set this to the current sea level pressure in the area for correct altitude readings
-    Pressure defaultSeaLevelPressure = WeatherHelper.MeanSeaLevel;
+    Pressure defaultSeaLevelPressure = WeatherHelper.MeanSeaLevelPressure;
     Length stationHeight = Length.FromMeters(640); // Elevation of the sensor
     var ftI2cBus = ft232h.CreateOrGetI2cBus(ft232h.GetDefaultI2cBusNumber());
     var i2cDevice = ftI2cBus.CreateDevice(Bmp280.DefaultI2cAddress);