極座標下的histogram2d

怎麼用PYTHON畫這種圖

主要代碼如下，最後附有我的原始代碼。下面這個圖可以用來輔助看

__author__ = 'kong'import matplotlib.pyplot as pltimport matplotlib.tri as triimport numpy as npimport mathimport netCDF4 as ncimport globimport mathimport matplotlibfrom mpl_toolkits.basemap import Basemap, shiftgrid#,cmfrom matplotlib import cmimport matplotlib.pyplot as pltfrom matplotlib.colors import LogNormimport netCDF4 as ncimport numpy as npfrom matplotlib import gridspecimport mathimport matplotlib.tri as triimport matplotlib.colors as matcolorprint("import finish")##-- circle and angle size --#### max_r min_r是畫圖時風速顯示的最大最小值## n 代表將風速大小，風向分為多少份，當然他們的份數可以不同，不過在這裏，它們的份數是一樣的## nn是n-1，沒甚麼意義，只是方便某些地方而已n=30#2*math.pi/0.01 max_r = 50max_r=float(max_r)min_r = 0min_r=float(min_r)nn=n-1##--read data--##ncfile=nc.Dataset(ifiles)u=ncfile.variables["u"][:]v=ncfile.variables["v"][:]ncfile.close()##--analyse data--##uu=u[np.where(u<10000)]vv=v[np.where(v<10000)]wind=np.sqrt(uu**2+vv**2)d=np.angle(uu+vv*1j)##--get draw data--####-- a is nn*nn, b and c is n --####b是風速大小的分界點##c是角度大小的分界點##詳情可以查看histogram2d函數a,b,c=np.histogram2d(wind,d,[np.arange(min_r,max_r,(max_r-min_r)/n),np.arange(-math.pi,math.pi,2*math.pi/n)])##將a一維化，它本來是二維的大小是(nn,nn)，一維化後大小是nnaa=a.flatten()##為了使aa和bbcc大小一樣，需要對它們做處理bb=(b[0:-1]+b[1:])/2cc=(c[0:-1]+c[1:])/2bbb=np.repeat(bb[np.newaxis,...],nn)ccc=np.zeros(nn*nn)for i in range(nn):  ccc[nn*i:nn*i+nn]=cc[:]##--start plot--##np.set_printoptions(threshold='nan')##--投影到極座標系下---##triang = tri.Triangulation(bbb*np.cos(ccc),bbb*np.sin(ccc))plt.tricontourf(triang, aa,levels=np.arange(0,np.max(aa)*1.1,1000))plt.colorbar()##--draw Auxiliary line--##my_r=np.arange(0,50,10)for i in my_r:  plt.plot(i*np.cos(ccc[0:n]),i*np.sin(ccc[0:n]),'r-')plt.plot([-70,70],[0,0],'r')plt.plot([0,0],[-70,70],'r')##--draw average (the white line)--##avg_wind=np.zeros(np.shape(cc))i_cc=0for i_angle in cc:  avg_wind[i_cc]=np.sum(a[:,i_cc]*bb)/np.sum(a[:,i_cc])  i_cc=i_cc+1new_avg_wind=np.zeros(np.shape(cc)[0]+1)new_cc=np.zeros(np.shape(cc)[0]+1)new_avg_wind[0:-1]=avg_windnew_cc[0:-1]=ccnew_avg_wind[-1]=avg_wind[0]new_cc[-1]=cc[0]plt.plot(new_avg_wind*np.cos(new_cc),new_avg_wind*np.sin(new_cc),'w-')plt.title("TITLE")plt.gca().set_aspect('equal')plt.savefig("meow.png", dpi=1000)#max dpi2300

下面是我的原代碼

__author__ = 'kong'import matplotlib.pyplot as pltimport matplotlib.tri as triimport numpy as npimport mathimport netCDF4 as ncimport globimport mathprint("import finish")n=30#2*math.pi/0.01max_r = 50max_r=float(max_r)min_r = 0min_r=float(min_r)nn=n-1outputa=np.zeros((n-1,n-1))files=glob.glob("//home//konghoiio//mydata//khidata//smallarea//*//my1area*.nc")numfiles=np.shape(files)i=0for ifiles in files:  i=i+1  print("i=%d, numfiles=%d"%(i,numfiles[0]))  b=nc.Dataset(ifiles)  u=b.variables["u"][:]  v=b.variables["v"][:]  b.close()  uu=u[np.where(u<10000)]  vv=v[np.where(v<10000)]  wind=np.sqrt(uu**2+vv**2)  d=np.angle(uu+vv*1j)  rose_aa,rose_bb,rose_cc=np.histogram2d(wind,d,[np.arange(min_r,max_r,(max_r-min_r)/n),np.arange(-math.pi,math.pi,2*math.pi/n)])  outputa=outputa+rose_aa  del(rose_aa)  del(d)  del(wind)  del(uu)  del(vv)  del(u)  del(v)da=nc.Dataset("windrosedata","w",format="NETCDF4")da.createDimension("n",n)da.createDimension("nn",nn)da.createVariable("a","f8",("nn","nn"))da.createVariable("b","f8",("n",))da.createVariable("c","f8",("n",))da.variables["a"][:]=outputada.variables["b"][:]=rose_bbda.variables["c"][:]=rose_ccda.description="wind rose data"da.author="konghoiio"da.createdate="2016-09-26"da.data_maxr=max_rda.data_minr=min_rda.data_n=nda.data_nn=nnda.close()print("FINISH")

上面的處理數據，下面是畫圖

__author__ = 'kong'import matplotlibfrom mpl_toolkits.basemap import Basemap, shiftgrid#,cmfrom matplotlib import cmimport matplotlib.pyplot as pltfrom matplotlib.colors import LogNormimport netCDF4 as ncimport numpy as npfrom matplotlib import gridspecimport mathimport matplotlib.tri as triimport matplotlib.colors as matcolordef drawwindrose():  filename="data//lljwindrosedatafor567"  plttitle="950hPa Wind Rose when LLJ for MJJ"  da=nc.Dataset(filename)  a=da.variables["a"][:]  b=da.variables["b"][:]  c=da.variables["c"][:]  n=da.data_n  nn=da.data_nn  da.close()  '''  da=nc.Dataset("data//lljwindrosedata")  aaa=da.variables["a"][:]  da.close()  a=a0-aaa  '''  aa=a.flatten()  bb=(b[0:-1]+b[1:])/2  cc=(c[0:-1]+c[1:])/2  bbb=np.repeat(bb[np.newaxis,...],nn)  ccc=np.zeros(nn*nn)  for i in range(nn):    ccc[nn*i:nn*i+nn]=cc[:]  np.set_printoptions(threshold='nan')  triang = tri.Triangulation(bbb*np.cos(ccc),bbb*np.sin(ccc))  plt.tricontourf(triang, aa,levels=np.arange(0,np.max(aa)*1.1,1000))  plt.colorbar()  my_r=np.arange(0,50,10)  for i in my_r:    plt.plot(i*np.cos(ccc[0:n]),i*np.sin(ccc[0:n]),'r-')  plt.plot([-70,70],[0,0],'r')  plt.plot([0,0],[-70,70],'r')  avg_wind=np.zeros(np.shape(cc))  i_cc=0  for i_angle in cc:    avg_wind[i_cc]=np.sum(a[:,i_cc]*bb)/np.sum(a[:,i_cc])    i_cc=i_cc+1  new_avg_wind=np.zeros(np.shape(cc)[0]+1)  new_cc=np.zeros(np.shape(cc)[0]+1)  new_avg_wind[0:-1]=avg_wind  new_cc[0:-1]=cc  new_avg_wind[-1]=avg_wind[0]  new_cc[-1]=cc[0]  plt.plot(new_avg_wind*np.cos(new_cc),new_avg_wind*np.sin(new_cc),'w-')  plt.title(plttitle)  plt.gca().set_aspect('equal')  #plt.savefig("20160921//"+filename[6:-4]+".png", dpi=1000)#max dpi2300  plt.savefig("20160921//"+plttitle+".png", dpi=1000)#max dpi2300

                                             
        0        0           
	
				
	
					   
極座標下的histogram2d
	  	   
的
	  	   
的
	  	   
的
	  	   
的
	  	   
的
	  	   
的
	  	   
的
	  	   
的
	  	   
的
	  	   
的
	  	   
的
	  	   
的
	  	   
的
	  	   
的
	  	   
的
	  	   
的
	  	   
的
	     		  
	  	   
浅谈StringBuffer类初始容量及扩容
	  	   
11E求最大公约数和最小公倍数
	  	   
Hadoop性能调优
	  	   
JRE、JDK和SDK分别是什么
	  	   
Visual Studio (VS) 操作小技巧（持续更新）
	  	   
極座標下的histogram2d
	  	   
腾讯云服务器上搭建mysql 如何使云数据库能外网访问
	  	   
Java 自定义标签
	  	   
实战 ：Spring MVC + 注解 +SqlServer 框架搭建及详解
	  	   
黑盒测试和白盒测试
	  	   
ios获取所有相册的视频并播放
	  	   
XML文件
	  	   
codeforces 733D. Kostya the Sculptor
	  	   
akka打包可执行jar错误