Code
Module
import sys
import numpy as N
import shapefile
import scipy.stats as stats
import pylab as P
from matplotlib.collections import PatchCollection
from pylab import Polygon
from mpl_toolkits.basemap import Basemap
P.style.use('bmh')
params = {'legend.fontsize': 8,'font.family': 'serif'}
P.rcParams.update(params)
Reading
tmp=N.genfromtxt('../shp/GP.csv',names=True,comments='#',delimiter=';',dtype=None,encoding='utf-8')
pds = tmp['Name']
file = '../dat/gwlneudatum.dat'
dd=N.genfromtxt(file,usecols=(0),skip_header=0,dtype="I")
mm=N.genfromtxt(file,usecols=(1),skip_header=0,dtype="I")
jj=N.genfromtxt(file,usecols=(2),skip_header=0,dtype="I")
gw=N.genfromtxt(file,usecols=(3),skip_header=0,dtype="S")
id = N.where((jj>=2015)&(jj<=2019))[0]
dd = dd[id]
mm = mm[id]
jj = jj[id]
gw = gw[id]
go = N.array(list(set(gw)))
mo = N.arange(1,13,1)
go = N.array(['WA','HNFZ','NWZ'])#N.array(['WZ','BM','TRM','HM','SWZ','TRW','TM','WA','HNFZ','NWZ'])
print (go)
print (mo)
Plotting
filter = 'Filter1'
paras = {'bahn':[0,30,'RdPu','[N]'],'tmax':[0,30,'Spectral_r',r'[$^\circ$C]'],'nied':[0,8,'YlGnBu','[mm/d]'],'wmax':[5,20,'BuGn','[m/s]']}
for gx in go:
for mx in mo:
jd = N.where((mm==mx)&(gw==gx))[0]
if(len(jd)>0.):
print (gx,mx,len(jd))
fig = P.figure(figsize=(12,4))
i = 0
#for para in paras:
for para in ['bahn','tmax','nied','wmax']:
i = i+1
ax = fig.add_subplot(1,4,i)
m = Basemap(projection='merc',llcrnrlat=47,urcrnrlat=55.2,llcrnrlon=5,urcrnrlon=16,resolution='l')
m.readshapefile('../shp/OSS_M1_PROD_DURCH_FL_EXT_geo', 'pds',zorder=10)
cmap=P.get_cmap(paras[para][2],15)
norm=P.Normalize(paras[para][0],paras[para][1])
patches = []
avg = []
for info, shape in zip(m.pds_info, m.pds):
pd = info['NAME']
if(para!='bahn'):
file = '../csv/obs-dwd/'+pd+'.csv'
data=N.genfromtxt(file,names=True,comments='#',delimiter=';',dtype=None,encoding='utf-8')
id = N.where((data['ja']>=2015)&(data['ja']<=2019))[0]
dum = data[para][id]
tmp = N.mean(dum[jd])
avg.append(tmp)
else:
file = '../csv/pds-uas/'+pd+'.csv'
data=N.genfromtxt(file,names=True,comments='#',delimiter=';',dtype=None,encoding='utf-8')
id = N.where((data['ja']>=2015)&(data['ja']<=2019))[0]
dum = data[filter][id]
tmp = N.sum(dum[jd])/len(jd)
#dum = []
#for ua in ['UA32','UA04','UA57','UA44','UA38','UA35','UA29']:
# xxx = data[ua][id]
# dum.append(xxx[jd])
#dum = N.array(dum)
#tmp = N.sum(N.sum(dum))/len(jd)
avg.append(tmp)
color=cmap(norm(tmp))
patches.append( Polygon(N.array(shape), True, color=color) )
pc = PatchCollection(patches, match_original=True, edgecolor='k', linewidths=0.1, zorder=2)#, alpha=0.5)
ax.add_collection(pc)
ax.axis('off')
Nshp=10
CM = P.scatter(N.zeros(Nshp),N.zeros(Nshp),s=N.zeros(Nshp),c=N.zeros(Nshp),cmap=cmap,vmin=paras[para][0],vmax=paras[para][1])
CB = m.colorbar(CM,location='bottom',extend='both',shrink=0.9,pad=0.05)
if(para=='bahn'):
CB.set_label('Meldungen '+paras[para][3],fontsize=12,weight='bold')
else:
CB.set_label(para+' '+paras[para][3],fontsize=12,weight='bold')
CB.ax.tick_params(labelsize=8,labelrotation=45)#,ha="center")#CB.ax.set_xticklabels(fontsize=10,rotation =45,ha="center")
avg = N.array(avg)
x,y = m(6,54.5)
P.text(x,y,r'$\varnothing$ %.1f'%N.mean(avg),fontsize=10,weight='bold')
#P.title('2015-2019: %02i-%s (N=%i)'%(mx,gx,len(jd)),fontsize=10,weight='bold')
fig.suptitle('Kritikalitaet von Wettermustern 2015-2019: %02i-%s (N=%i,%s)'%(mx,gx,len(jd),filter),fontsize=16,weight='bold',y=1.05)
fig.tight_layout()
fig.savefig('./img/%02i-%s_%s.png'%(mx,gx,filter),dpi=240,transparent=False,bbox_inches='tight',pad_inches=0.0)