Verify GFS dewpoint with MesoWest observation for class assignment

verify_GFS_pygrib.py

# Brian Blaylock
# ATMOS 6910 Homework
# 10 November 2016

# Look at the figures here: 
# http://kbkb-wx-python.blogspot.com/2016/11/verifying-gfs-dewpoint-data-with.html

"""
[X] Reads a GRiB forecast file from the GFS (15%)
[X] Extracts Temperature and Relative Humidity from ingest (5%)
[X] Calculates dew point from temperature and RH (10%)
[X] Plots All three maps (15%)
[X] Use the Mesowest API to verify the dew point programmatically (15%)
[X] Create unit tests to verify your dew point function (10%)
[X] Use all the best practices (Use the ASTG as your guide.) (15%)
[X] Include quality documentation (15%)
"""

import numpy as np
import pygrib
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

import sys
sys.path.append('/uufs/chpc.utah.edu/common/home/u0553130/pyBKB_v2/')
sys.path.append('B:/pyBKB_v2/')
from BB_MesoWest.MesoWest_radius import get_mesowest_radius


def K_to_C(K):
    """
    Convert a temperature in Kelvin to Celsius
    """
    return K-273.15


def dwptRH_to_Temp(dwpt, RH):
    """
    Convert a dew point temerature and relative humidity to an air temperature.
    Equation from:
    http://andrew.rsmas.miami.edu/bmcnoldy/humidity_conversions.pdf

    Input:
        dwpt - Dew point temperature in Celsius
        RH - relative humidity in %
    Output:
        Temp - Temperature in Celsius
    """
    a = 17.625
    b = 243.04
    Temp = b * (a*dwpt/(b+dwpt)-np.log(RH/100.)) / (a+np.log(RH/100.)-(a*dwpt/(b+dwpt)))
    return Temp


def Tempdwpt_to_RH(Temp, dwpt):
    """
    Convert a temperature and dew point temerature to relative humidity.
    Equation from:
    http://andrew.rsmas.miami.edu/bmcnoldy/humidity_conversions.pdf

    Input:
        Temp - Temperature in Celsius
        dwpt - Dew point temperature in Celsius

    Output:
        RH - relative humidity in % 
    """
    a = 17.625
    b = 243.04
    RH = 100*(np.exp((a*dwpt/(b+dwpt)))/np.exp((a*Temp/(b+Temp))))
    return RH


def TempRH_to_dwpt(Temp, RH):
    """
    Convert a temperature and relative humidity to a dew point temperature.
    Equation from:
    http://andrew.rsmas.miami.edu/bmcnoldy/humidity_conversions.pdf

    Input:
        Temp - Air temperature in Celsius
        RH - relative humidity in %
    Output:
        dwpt - Dew point temperature in Celsius 
    """
    # Check if the Temp coming in is in celsius and if RH is between 0-100%
    passed = False
    test_temp = temp<65
    
    if np.sum(test_temp) == np.size(temp):
        passed = True
        test_rh = np.logical_and(RH<=100,RH>=0)
        if np.sum(test_rh) == np.size(RH):
            passed = True
        else:
            print "faied relative humidity check"
    else:
        print "faild temperature check"

    if passed == True:
        a = 17.625
        b = 243.04
        dwpt = b * (np.log(RH/100.) + (a*Temp/(b+Temp))) / (a-np.log(RH/100.)-((a*Temp)/(b+Temp)))
        return dwpt

    else:
        print "TempRH_to_dwpt input requires a valid temperature and humidity."
        return "Input needs a valid temperature (C) and humidity (%)."


def draw_world_map():
    """
    Draw a custom basemap
    """
    ## Draw Background basemap
    bot_left_lat = -90
    bot_left_lon = 0
    top_right_lat = 90
    top_right_lon = 360

    ## Map in cylindrical projection (data points may apear skewed)
    m = Basemap(resolution='l', projection='cyl', \
        llcrnrlon=bot_left_lon, llcrnrlat=bot_left_lat, \
        urcrnrlon=top_right_lon, urcrnrlat=top_right_lat)
    return m


def draw_utah_map():
    """
    Draw a custom basemap
    """
    ## Draw Background basemap
    bot_left_lat  =36.5
    bot_left_lon  =-114.5
    top_right_lat =42.5
    top_right_lon = -108.5

    ## Map in cylindrical projection (data points may apear skewed)
    m = Basemap(resolution='l', projection='cyl', \
        llcrnrlon=bot_left_lon, llcrnrlat=bot_left_lat, \
        urcrnrlon=top_right_lon, urcrnrlat=top_right_lat)
    return m


def pluck_point(stn_lat,stn_lon,grid_lat,grid_lon):
    """
    Get the WRF data for the point nearest the MesoWest station
    Figure out the nearest lat/lon in the HRRR domain for the station location
    Input:
        stn_lat, stn_lon - The latitude and longitude of the station
        grid_lat, grid_lon - The 2D arrays for both latitude and longitude
    Output:
        The index of the flattened array
    """

    abslat = np.abs(grid_lat-stn_lat)
    abslon = np.abs(grid_lon-stn_lon)
    # Element-wise maxima. Plot this with pcolormesh to see what I've done.
    c = np.maximum(abslon, abslat)
    # The minimum maxima (which which is the nearest lat/lon
    latlon_idx = np.argmin(c)
    return latlon_idx


if __name__ == '__main__':
    # Open the grib file and grab the variables we need.
    BASE = '/uufs/chpc.utah.edu/common/home/u0079358/public_html/atmos6910/Class_4/'
    FILE = 'fh.0003_tl.press_gr.0p50deg'
    grbs = pygrib.open(BASE + FILE)

    TEMP = grbs.select(name='2 metre temperature')[0]
    RH = grbs.select(name='Relative humidity')[31]
    print 'Grabbed:', TEMP
    print 'Grabbed:', RH
    temp, lat, lon = TEMP.data()
    rh, lat, lon = RH.data()
    DATE = TEMP.validDate

    # Convert temperature to Celsius and then find the dew point.
    temp = K_to_C(temp)
    dwpt = TempRH_to_dwpt(temp, rh)
    # For Utah Maps, need to fix the lon for West longitude values
    lon2 = np.copy(lon)
    lon2[lon > 180] = lon[lon > 180] - 360

    # Compare the GFS model with MesoWest observations.
    # 1) Get the mesowest data for a radius around KSLC for the same time assert
    #    the grib2 file.
    a = get_mesowest_radius(DATE, '10', extra='&state=ut')

    # 2) For each station, pluck the corresponding point from the GFS grid.
    point_dwpt = np.array([])
    point_lat = np.array([])
    point_lon = np.array([])
    for i in range(0, len(a['STID'])):
        # Plucked index (of a flattened array)
        pi = pluck_point(a['LAT'][i], a['LON'][i], lat, lon2)
        # Append with the plucked Value
        point_dwpt = np.append(point_dwpt, dwpt.flat[pi])
        point_lat = np.append(point_lat, lat.flat[pi])
        point_lon = np.append(point_lon, lon2.flat[pi])

    # 3) Finally, calcualte the difference between the GFS and MesoWest
    #    dewpoints
    dwpt_diff = point_dwpt-a['dew_point_temperature']


    # --- Make some totally awesome figures!-----------------------------------
    # Plot world maps of the data.
    m = draw_world_map()
    mUtah = draw_utah_map()

    # World map of surface temperature
    plt.figure(1)
    m.pcolormesh(lon, lat, temp, cmap='afmhot_r', vmin=-20, vmax=40)
    m.drawcoastlines()
    m.drawcountries()
    m.drawstates(linewidth=.25)
    plt.colorbar(shrink=.5)
    plt.title('2-m Temperature (C)')
    plt.savefig('world_temp.png', bbox_inches='tight')

    # World map of surface relative humidity
    plt.figure(2)
    m.drawcoastlines()
    m.drawcountries()
    m.drawstates(linewidth=.25)
    m.pcolormesh(lon, lat, rh, cmap='BrBG', vmin=0, vmax=100)
    plt.colorbar(shrink=0.5)
    plt.title('2-m Relative Humidity (%)')
    plt.savefig('world_rh.png', bbox_inches='tight')

    # World map of surface dew point temperature
    plt.figure(3)
    m.drawcoastlines()
    m.drawcountries()
    m.drawstates(linewidth=.25)
    m.pcolormesh(lon, lat, dwpt, cmap='BrBG', vmin=-20, vmax=40)
    plt.colorbar(shrink=.5)
    plt.title('2-m Dew Point Temperature (C)')
    plt.savefig('world_dwpt.png', bbox_inches='tight')

    # Utah map of GFS and MesoWest dewpoints
    plt.figure(4)
    mUtah.drawcoastlines()
    mUtah.drawcountries()
    mUtah.drawstates(linewidth=.25)
    mUtah.pcolormesh(lon2, lat, dwpt, cmap='BrBG', vmin=-20, vmax=40)
    mUtah.scatter(a['LON'], a['LAT'], c=a['dew_point_temperature'], cmap='BrBG', vmax=40, vmin=-20)
    plt.colorbar(shrink=.5)
    plt.title('Dew Point Temperature (C)')
    plt.savefig('utah_dwpt.png', bbox_inches='tight')

    # Utah map of difference between GFS and Mesowest dewpoints (GFS-MesoWest)
    plt.figure(5)
    mUtah.drawcoastlines()
    mUtah.drawcountries()
    mUtah.drawstates(linewidth=.25)
    mUtah.scatter(a['LON'], a['LAT'], c=dwpt_diff, cmap='bwr', vmax=3, vmin=-3)
    plt.colorbar(shrink=.5)
    plt.title('Dew Point Differece (C): GFS - Station')
    plt.savefig('utah_dwpt_diff.png', bbox_inches='tight')

    # Utah map of MesoWest stations and GFS data point used for the comparison.
    plt.figure(6)
    mUtah.drawcoastlines()
    mUtah.drawcountries()
    mUtah.drawstates(linewidth=.25)
    plt.title('Station Verification Point')
    # Draw line between each station and the verification point
    for i in range(0, len(a['STID'])):
        mUtah.plot([a['LON'][i], point_lon[i]], [a['LAT'][i], point_lat[i]], color='k')
    # verification grid point and stations on a map
    mUtah.scatter(point_lon, point_lat, s=70, c='blue')
    mUtah.scatter(a['LON'], a['LAT'], s=40, c='red')
    plt.savefig('utah_ver_points.png', bbox_inches='tight')

    plt.show(block=False)