#!/usr/bin/python2 # # (c) author Syn Gron (syngron(a)gmail.com) # # Data provided by German Weather Service (Deutscher Wetterdienst, DWD): http://www.dwd.de # Monthly data available from: http://www.dwd.de/bvbw/appmanager/bvbw/dwdwwwDesktop?_nfpb=true&_pageLabel=_dwdwww_klima_umwelt_klimadaten_deutschland&T82002gsbDocumentPath=Navigation%2FOeffentlichkeit%2FKlima__Umwelt%2FKlimadaten%2Fkldaten__kostenfrei%2Fkldat__D__mittelwerte__node.html%3F__nnn%3Dtrue # Sonnenscheindauer 1981 - 2010 (aktueller Standort) from mpl_toolkits.basemap import Basemap, cm import numpy as np import matplotlib.pyplot as plt import csv import matplotlib import scipy.interpolate import sys import calendar import locale #np.set_printoptions(threshold=5) # minimal and maximal days for common scaling data_min_fix=0 data_max_fix=300 # resolution of the grid grd_res=0.01 #grd_res=1 # resolution of the map map_res='h' #map_res='c' # individual scaling or common scaling of the months individual=True #individual=False globalmin=sys.float_info.max globalmax=0 def setMap(lats, lons, data, showstations): x = np.array(lons) y = np.array(lats) z = np.array(data) # corners of the map lon_min = np.min(x)-1 lon_max = np.max(x)+1 lat_min = np.min(y)-1 lat_max = np.max(y)+1 if individual: data_min = np.min(z) data_max = np.max(z) else: data_min = data_min_fix data_max = data_max_fix # general map m = Basemap(projection = 'merc',llcrnrlat=lat_min, urcrnrlat=lat_max,llcrnrlon=lon_min, urcrnrlon=lon_max,resolution=map_res, area_thresh=100) m.drawcoastlines() m.drawstates() m.drawcountries() if not showstations: y_inc = (lat_max - lat_min) / grd_res x_inc = (lon_max - lon_min) / grd_res y_steps = np.linspace(lat_min, lat_max + grd_res, y_inc) x_steps = np.linspace(lon_min, lon_max + grd_res, x_inc) x_steps, y_steps = np.meshgrid(x_steps, y_steps) zgrd = matplotlib.mlab.griddata(x, y, z, x_steps, y_steps) xgrd, ygrd = m.makegrid(zgrd.shape[1], zgrd.shape[0]) #xgrd,ygrd= np.meshgrid(x_steps,y_steps) #zgrd = scipy.interpolate.griddata((x, y), z, (xgrd, ygrd), method='nearest') xgrd,ygrd = m(xgrd,ygrd) m.contourf(xgrd,ygrd,zgrd,cmap=plt.cm.hot_r,vmin=data_min, vmax=data_max) else: # measurement stations xpts,ypts = m(lons,lats) m.plot(xpts,ypts,'bo') #for xpt, ypt, caption in zip(xpts, ypts, captions): # plt.text(xpt,ypt,caption, fontsize=5) fig, axes = plt.subplots(figsize=(18,10)) # otherwise month names will be in locale language locale.setlocale(locale.LC_ALL, 'en_US') myrange=13 lats=[] lons=[] captions=[] alldata=[[] for i in range(myrange)] with open('DWD_Sonnenschein.csv', 'rb') as f: reader = csv.reader(f) rownum=0 for row in reader: if rownum>2: lats.append(float(row[4][:2]) + float(row[4][5:7])/60.0) lons.append(float(row[5][:2]) + float(row[5][5:7])/60.0) captions.append(str(row[0])) for i in range(myrange): alldata[i].append(float(row[i+7])) rownum+=1 for i in range(myrange): print i data=alldata[i] if i == 12: data = np.divide(data,12) # yearly data as mean tmpmax=np.max(data) tmpmin=np.min(data) if tmpmax > globalmax: globalmax=tmpmax if tmpmin < globalmin: globalmin=tmpmin plt.subplot(3,6,i+1) if i == 12: plt.title('Avg. per year') else: plt.title(calendar.month_name[i+1]) setMap(lats, lons, data, False) if individual: plt.colorbar() plt.subplot(3,6,14) plt.title('Measurement stations') setMap(lats, lons, data, True) print globalmax print globalmin if not individual: #ax, _ = matplotlib.colorbar.make_axes(plt.gca()) cax = fig.add_axes([0.92, 0.1, 0.01, 0.8]) cbar = matplotlib.colorbar.ColorbarBase(cax, cmap=plt.cm.hot_r, norm=matplotlib.colors.Normalize(vmin=data_min_fix, vmax=data_max_fix)) cbar.set_label('Hours of sunlight') cbar.set_clim(data_min_fix, data_max_fix) fig.canvas.set_window_title('Monthly hours of sunlight') #if individual: # plt.savefig('hours_of_sunlight_individual_scale.png', dpi=300) #else: # plt.savefig('hours_of_sunlight.png', dpi=300) plt.show()