Unidata / MetPy

MetPy is a collection of tools in Python for reading, visualizing and performing calculations with weather data.
https://unidata.github.io/MetPy/
BSD 3-Clause "New" or "Revised" License
1.26k stars 415 forks source link

Calculating specific humidity from dew point temperature #791

Closed avatar101 closed 5 years ago

avatar101 commented 6 years ago

Hi,

I would like to request a function to calculate specific humidity directly from dew point temperature. Also, if the function considers the option to use Latent Heat for Ice to make sure to calculate specific humidity correctly in the polar region. I could share the code of the function I wrote but, it doesn't make use of pint quantities.

Thanks!

jrleeman commented 6 years ago

Sharing is great! We're happy to help adapt code or do it ourselves when we get to this!

dopplershift commented 6 years ago

Certainly it's always helpful to start from at least a sketch that someone else has found useful--that might answer some design questions or save us from researching something ourselves.

avatar101 commented 6 years ago

Here is the python function I wrote to calculate specific humidity. The formula I'm using is referenced in the comments. In case of any further questions, please feel free to ask. Cheers,

import numpy as np
class sounding_df:

    def __init__(self, df):
        # this __init__ method is used to create new objects inside the class
        self.p = np.array(df['Pressure'])
        self.temp = np.array(df['Temp']) # in Kelvin
        self.dew = np.array(df['Dewpt']) # in Kelvin

    def spec_humidity(self, latent='water'):
        """Calculates SH automatically from the dewpt. Returns in g/kg"""
        # Declaring constants
        self.e0 = 6.113 # saturation vapor pressure in hPa
        # e0 and Pressure have to be in same units
        self.c_water = 5423 # L/R for water in Kelvin
        self.c_ice = 6139 # L/R for ice in Kelvin
        self.T0 = 273.15 # Kelvin

        if latent == 'water' or latent == 'Water':
            self.c = self.c_water    # using c for water

        else:
            self.c= self.c_ice       # using c_ice for ice, clear state

        # saturation vapor not required, uncomment to calculate it (units in hPa becuase of e0)
        #sat_vapor = self.e0 * np.exp((self.c * (self.temp -self.T0))/(self.temp * self.T0)) 

        #calculating specific humidity, q directly from dew point temperature
        #using equation 4.24, Pg 96 Practical Meteorolgy (Roland Stull)
        self.q = (622 * self.e0 * np.exp(self.c * (self.dew - self.T0)/(self.dew * self.T0)))/self.p # g/kg 
        # 622 is the ratio of Rd/Rv in g/kg
        return self.q
ytian0704 commented 3 years ago

I think there's one mistake in your function. Checking the equation 4.24 in the book, p should be in the front like (622*e0/p) as a whole instead of being outside of the brackets. I think the formula should be:

self.q = (622 self.e0/self.p) np.exp(self.c (self.dew - self.T0)/(self.dew self.T0))

where q in the unit of g/kg, e0 and p in the same unit in order to be cancelled out, c, dew and T0 in the unit of K

dcamron commented 3 years ago

@ytian0704 this was implemented differently in MetPy in the referenced PR #1093. If your correction is still helpful for @avatar101, then thanks for that!