DARL = 9.8C/1000m
In [1]:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import json
import gspread
from oauth2client.service_account import ServiceAccountCredentials
import pytz
from datetime import datetime
# df = pd.read_csv('sounding.csv', names= cnames, header=0)
def computeForecast(df):
# Altitude in feet
df['ALT'] = df.HEIGHT*3.28084
df.TEMP = df.TEMP/10
df.DEWPT = df.DEWPT/10
df['SPREAD'] = df.TEMP - df.DEWPT
df['VIRTT'] = ((df.TEMP + 273.15)/(1 - 0.378*(pow(6.11*10,7.5*df.DEWPT/(237.7+df.DEWPT))/df.PRESSURE)))-273.15
# DALR, T.I.
max_t = df.iloc[0,3]
surface_alt = df.iloc[0,2]
df['DALR'] = max_t - (9.8/1000)*(df.HEIGHT-surface_alt)
df['TI'] = df.VIRTT - df.DALR
# Scorer
df['TK'] = df.TEMP+273.15
df['Y'] = (df.TEMP-df.TEMP.shift(-1))/(df.HEIGHT.shift(-1)-df.HEIGHT)
df['WV_MPS'] = df.WIND_SPD*0.514
df['L2'] = (((0.00986-df.Y)/df.TK)*(9.81/pow(df.WV_MPS,2))-(1/4)*pow(((9.81/287-df.Y)/df.TK),2))*100000
return df
def processSummary(col, sdfs, hour):
dfs = pd.DataFrame()
for sdf in sdfs:
dfs['ALT'] = sdf['ALT']
dfs['UTC_' + str(hour)] = sdf[col]
hour = hour + 1
dfs = dfs.sort_values(by=['ALT'], ascending=[0])
return dfs
def numberToLetters(q):
q = q - 1
result = ''
while q >= 0:
remain = q % 26
result = chr(remain+65) + result;
q = q//26 - 1
return result
def uploadToGoogle(df, sheet):
# Upload the results to Google Drive
scope = ['https://spreadsheets.google.com/feeds']
credentials = ServiceAccountCredentials.from_json_keyfile_name('flightbit-key.json', scope)
gc = gspread.authorize(credentials)
wb = gc.open_by_key('1q8YhqClfmTF_be4QnqgDGMqJfDj_ELdQSrrjlqPb770')
ws = wb.worksheet(sheet)
# columns names
columns = df.columns.values.tolist()
# selection of the range that will be updated
cell_list = ws.range('A1:'+numberToLetters(len(columns))+'1')
# modifying the values in the range
for cell in cell_list:
val = columns[cell.col-1]
if type(val) is str:
val = val.decode('utf-8')
cell.value = val
# update in batch
ws.update_cells(cell_list)
# number of lines and columns
num_lines, num_columns = df.shape
# selection of the range that will be updated
cell_list = ws.range('A2:'+numberToLetters(num_columns)+str(num_lines+1))
# modifying the values in the range
for cell in cell_list:
val = df.iloc[cell.row-2,cell.col-1]
if type(val) is str:
val = val.decode('utf-8')
elif isinstance(val, (int, long, float, complex)):
# note that we round all numbers
val = val
cell.value = val
# update in batch
ws.update_cells(cell_list)
cnames = ["TYPE","PRESSURE","HEIGHT","TEMP","DEWPT","WIND_DIR","WIND_SPD"]
now = datetime.now(tz=pytz.utc)
df = pd.read_csv('http://rucsoundings.noaa.gov/get_soundings.cgi?data_source=Op40&latest=latest&n_hrs=24&' +
'fcst_len=shortest&airport=kawo&text=Ascii%20text%20%28GSD%20format%29&hydrometeors=false&start=latest',
names= cnames, skiprows=1, delim_whitespace=1)
hour = int(df.iloc[0,1])
print "Start Hour (UTC): " + str(hour)
df = df[df.TYPE.isin(['4', '5', '9'])]
df = df.reset_index(drop=True)
# Find the row indexes where new soundings start
idxs = df[df['TYPE'] == '9'].index.tolist()
df[cnames] = df[cnames].astype(float)
sdfs = []
for idx, val in enumerate(idxs):
if val != idxs[-1]:
start = val
end = idxs[idx+1]
sdfs.append(df[start:end].copy().reset_index(drop=True))
else:
sdfs.append(df[val:].copy().reset_index(drop=True))
# Process the soundings
for sdf in sdfs:
computeForecast(sdf)
# Generate summaries
for summary in ['TI', 'SPREAD', 'WIND_SPD', 'WIND_DIR', 'L2']:
df_s = processSummary(summary, sdfs, hour)
uploadToGoogle(df_s, summary)
In [9]:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import json
import gspread
from oauth2client.service_account import ServiceAccountCredentials
import pytz
from datetime import datetime
cnames = ["PRESSURE","HEIGHT","TEMP","DEWPT","WIND_DIR","WIND_SPD"]
df = pd.read_csv('sounding_kawo.csv', names= cnames, header=0)
def computeForecast(df):
# Altitude in feet
df['ALT'] = df.HEIGHT*3.28084
# Convert from T in 10th of degrees celcius
df.TEMP=df.TEMP
df.DEWPT=df.DEWPT
df['SPREAD'] = df.TEMP - df.DEWPT
df['VIRTT'] = ((df.TEMP + 273.15)/(1 - 0.378*(pow(6.11*10,7.5*df.DEWPT/(237.7+df.DEWPT))/df.PRESSURE)))-273.15
# DALR, T.I.
max_t = df.iloc[0,2]
surface_alt = df.iloc[0,1]
df['DALR'] = max_t - (9.8/1000)*(df.HEIGHT-surface_alt)
df['TI'] = df.VIRTT - df.DALR
# Scorer
df['TK'] = df.TEMP+273.15
df['Y'] = (df.TEMP-df.TEMP.shift(-1))/(df.HEIGHT.shift(-1)-df.HEIGHT)
df['WV_MPS'] = df.WIND_SPD*0.514
df['L2'] = (((0.00986-df.Y)/df.TK)*(9.81/pow(df.WV_MPS,2))-(1/4)*pow(((9.81/287-df.Y)/df.TK),2))*100000
return df
result = computeForecast(df)
result.to_csv('kawo_forecast.csv')
In [3]:
from mpl_toolkits.basemap import Basemap, cm
# requires netcdf4-python (netcdf4-python.googlecode.com)
from netCDF4 import Dataset as NetCDFFile
import numpy as np
import matplotlib.pyplot as plt
# plot rainfall from NWS using special precipitation
# colormap used by the NWS, and included in basemap.
nc = NetCDFFile('nws_precip_conus_20160116.nc')
# data from http://water.weather.gov/precip/
prcpvar = nc.variables['amountofprecip']
data = 0.01*prcpvar[:]
latcorners = nc.variables['lat'][:]
loncorners = -nc.variables['lon'][:]
lon_0 = -nc.variables['true_lon'].getValue()
lat_0 = nc.variables['true_lat'].getValue()
# create figure and axes instances
fig = plt.figure(figsize=(8,8))
ax = fig.add_axes([0.1,0.1,0.8,0.8])
# create polar stereographic Basemap instance.
m = Basemap(projection='stere',lon_0=lon_0,lat_0=90.,lat_ts=lat_0,\
llcrnrlat=latcorners[0],urcrnrlat=latcorners[2],\
llcrnrlon=loncorners[0],urcrnrlon=loncorners[2],\
rsphere=6371200.,resolution='l',area_thresh=10000)
# draw coastlines, state and country boundaries, edge of map.
m.drawcoastlines()
m.drawstates()
m.drawcountries()
# draw parallels.
parallels = np.arange(0.,90,10.)
m.drawparallels(parallels,labels=[1,0,0,0],fontsize=10)
# draw meridians
meridians = np.arange(180.,360.,10.)
m.drawmeridians(meridians,labels=[0,0,0,1],fontsize=10)
ny = data.shape[0]; nx = data.shape[1]
lons, lats = m.makegrid(nx, ny) # get lat/lons of ny by nx evenly space grid.
x, y = m(lons, lats) # compute map proj coordinates.
# draw filled contours.
clevs = [0,1,2.5,5,7.5,10,15,20,30,40,50,70,100,150,200,250,300,400,500,600,750]
cs = m.contourf(x,y,data,clevs,cmap=cm.s3pcpn)
# add colorbar.
cbar = m.colorbar(cs,location='bottom',pad="5%")
cbar.set_label('mm')
# add title
plt.title(prcpvar.long_name+' for period ending '+prcpvar.dateofdata)
plt.show()
In [ ]: