sbirch · February 11, 2015 08:03 · Feb 11, 2015
diff --git a/analysis.py b/analysis.py
@@ -0,0 +1,81 @@
+import pandas as pd
+import numpy as np
+import re
+import matplotlib.pyplot as plt
+from mpl_toolkits.basemap import Basemap
+
+def load_grid(path):
+    f = open(path, 'r')
+
+    headers = None
+    for line in f:
+        if line.startswith('Lat Lon'):
+            headers = re.split('\s+', line.strip())
+            break
+
+    records = []
+    for line in f:
+        records.append([float(x) for x in re.split('\s+', line.strip())])
+
+    df = pd.DataFrame(np.array(records), columns=headers)
+    assert len(df) == len(set(zip(df['Lat'], df['Lon'])))
+
+    return df
+
+
+# https://eosweb.larc.nasa.gov/cgi-bin/sse/global.cgi?email=skip@larc.nasa.gov
+
+# Insolation Incident On A Horizontal Surface (kWh/m^2/day)
+horizontal_radiation = load_grid('Horizontal_Radiation.txt')
+# Air Temperature at 10 m Above The Surface Of The Earth (deg C)
+air_temp = load_grid('AirTemp_22yr_T10M.txt')
+
+insolation_efficiency = []
+
+assert np.all(horizontal_radiation['Lon'] == air_temp['Lon'])
+assert np.all(horizontal_radiation['Lat'] == air_temp['Lat'])
+
+for i in xrange(len(horizontal_radiation)):
+    if horizontal_radiation.loc[i,'Lat'] >= 0:
+        stat = (air_temp.iloc[i,8]+273.0) / horizontal_radiation.iloc[i,8]
+    else:
+        stat = (air_temp.iloc[i,2]+273.0) / horizontal_radiation.iloc[i,2]
+
+    insolation_efficiency.append((
+        horizontal_radiation.loc[i,'Lat'],
+        horizontal_radiation.loc[i,'Lon'],
+        stat))
+
+insolation_efficiency = pd.DataFrame(np.array(insolation_efficiency), columns=['Lat', 'Lon', 'ieff'])
+
+lats = np.arange(-90,90,1)
+lons = np.arange(-180,180,1)
+lons, lats = np.meshgrid(lons,lats)
+
+fig = plt.figure()
+ax = fig.add_axes([0.05,0.05,0.9,0.9])
+
+m = Basemap(projection='kav7',lon_0=0,resolution='l')
+
+m.drawcoastlines(linewidth=0.25)
+m.drawcountries(linewidth=0.25)
+m.drawmapboundary(fill_color='#A0D9F1')
+
+data = insolation_efficiency['ieff'].values
+data = data.reshape((180, 360))
+
+im1 = m.pcolormesh(
+    lons, lats,
+    data,
+    shading='flat',
+    cmap=plt.cm.jet,
+    latlon=True)
+
+cb = m.colorbar(im1,"bottom", size="5%", pad="2%")
+
+# draw parallels and meridians, but don't bother labelling them.
+m.drawparallels(np.arange(-90.,99.,30.))
+m.drawmeridians(np.arange(-180.,180.,60.))
+
+ax.set_title('Air Temperature @10m ($K$) / Horizontal Insolation ($kWh/m^2/day$) in Summer Month')
+plt.show()
No results found