極座標下的histogram2d

怎麼用PYTHON畫這種圖

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

__author__ = 'kong'
import matplotlib.pyplot as plt
import matplotlib.tri as tri
import numpy as np
import math
import netCDF4 as nc
import glob
import math
import matplotlib
from mpl_toolkits.basemap import Basemap, shiftgrid#,cm
from matplotlib import cm
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
import netCDF4 as nc
import numpy as np
from matplotlib import gridspec
import math
import matplotlib.tri as tri
import matplotlib.colors as matcolor
print("import finish")
##-- circle and angle size --##
## max_r min_r是畫圖時風速顯示的最大最小值
## n 代表將風速大小，風向分為多少份，當然他們的份數可以不同，不過在這裏，它們的份數是一樣的
## nn是n-1，沒甚麼意義，只是方便某些地方而已
n=30#2*math.pi/0.01 
max_r = 50
max_r=float(max_r)
min_r = 0
min_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)，一維化後大小是nn
aa=a.flatten()
##為了使aa和bbcc大小一樣，需要對它們做處理
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[:]

##--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=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("TITLE")
plt.gca().set_aspect('equal')
plt.savefig("meow.png", dpi=1000)#max dpi2300

下面是我的原代碼

__author__ = 'kong'
import matplotlib.pyplot as plt
import matplotlib.tri as tri
import numpy as np
import math
import netCDF4 as nc
import glob
import math
print("import finish")


n=30#2*math.pi/0.01
max_r = 50
max_r=float(max_r)
min_r = 0
min_r=float(min_r)
nn=n-1

outputa=np.zeros((n-1,n-1))

files=glob.glob("//home//konghoiio//mydata//khidata//smallarea//*//my1area*.nc")
numfiles=np.shape(files)
i=0
for 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"][:]=outputa
da.variables["b"][:]=rose_bb
da.variables["c"][:]=rose_cc

da.description="wind rose data"
da.author="konghoiio"
da.createdate="2016-09-26"
da.data_maxr=max_r
da.data_minr=min_r
da.data_n=n
da.data_nn=nn


da.close()

print("FINISH")

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

__author__ = 'kong'
import matplotlib
from mpl_toolkits.basemap import Basemap, shiftgrid#,cm
from matplotlib import cm
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
import netCDF4 as nc
import numpy as np
from matplotlib import gridspec
import math
import matplotlib.tri as tri
import matplotlib.colors as matcolor


def 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