Here is the output figures and my code.
The output:
gist:
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| # 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) |
