基於python的dea繪圖（dea模型python）

• 1、python的作圖包有哪些
• 2、Python如何運用matplotlib庫繪製3D圖形
• 3、Python作圖程序
• 4、Python氣象數據處理與繪圖(18)：泰勒圖
• 5、如何使用Python繪製GWAS分析中的曼哈頓圖和QQ圖
• 6、Python氣象數據處理與繪圖(2)：常用數據計算方法

1. matplotlib

該python繪圖包與matlab的繪圖功能類似

2. seaborn

用來進行統計數據可視化的工具包，繪製的圖像非常漂亮，該python package也是基於matplotlib的，是純粹由python開發的

3D圖形在數據分析、數據建模、圖形和圖像處理等領域中都有著廣泛的應用，下面將給大家介紹一下如何在Python中使用 matplotlib進行3D圖形的繪製，包括3D散點、3D表面、3D輪廓、3D直線（曲線）以及3D文字等的繪製。

準備工作：

python中繪製3D圖形，依舊使用常用的繪圖模塊matplotlib，但需要安裝mpl_toolkits工具包，安裝方法如下：windows命令行進入到python安裝目錄下的Scripts文件夾下，執行： pip install –upgrade matplotlib即可；Linux環境下直接執行該命令。

安裝好這個模塊後，即可調用mpl_tookits下的mplot3d類進行3D圖形的繪製。

下面以實例進行說明。

1、3D表面形狀的繪製

這段代碼是繪製一個3D的橢球表面，結果如下：

2、3D直線（曲線）的繪製

這段代碼用於繪製一個螺旋狀3D曲線，結果如下：

3、繪製3D輪廓

繪製結果如下：

相關推薦：《Python視頻教程》

4、繪製3D直方圖

繪製結果如下：

5、繪製3D網狀線

繪製結果如下：

6、繪製3D三角面片圖

繪製結果如下：

7、繪製3D散點圖

繪製結果如下：

實戰小程序：畫出y=x^3的散點圖

樣例代碼如下：

[python] view plain copy

#coding=utf-8

import pylab as y    #引入pylab模塊

x = y.np.linspace(-10, 10, 100)  #設置x橫坐標範圍和點數

y.plot(x, x*x*x,’or’)  #生成圖像

ax = y.gca()

ax.spines[‘right’].set_color(‘none’)

ax.spines[‘top’].set_color(‘none’)

ax.xaxis.set_ticks_position(‘bottom’)

ax.spines[‘bottom’].set_position((‘data’, 0))

ax.yaxis.set_ticks_position(‘left’)

ax.spines[‘left’].set_position((‘data’, 0))

ax.set_yticks([-1000, -500, 500, 1000])

y.xlim(x.min() , x.max() ) #將橫坐標設置為x的最大值和最小值

y.show() #顯示圖像

[python] view plain copy

import pylab as y

程序中引入的pylab屬於matplotlib的一個模塊，將其名字用y代替，其中包括了許多NumPy和pyplot模塊中常用的函數，方便用戶快速進行計算和繪圖，十分適合在IPython互動式環境中使用。

[python] view plain copy

y.np.linspace(-10, 10, 100)

此為numpy中的一個函數，返回的是等間距的值，numpy.linspace(a,b,c):a指的是開始位置，b表示的是結束位置，c表示產生點的個數（默認為50）

舉例：

[python] view plain copy

 np.linspace(2.0, 3.0, num=5)

array([ 2.  ,  2.25,  2.5 ,  2.75,  3.  ])

[python] view plain copy

y.plot(x, x*x*x,’or’)  #生成圖像

後面加上『o’表示為散點圖

‘r’可設置顏色為紅色，基本上和matlab的操作很像。

[python] view plain copy

y.xlim(x.min(), x.max())

這條語句使用了xlim函數，將橫坐標設置為x的大小

泰勒圖繪製的核心思想是設計一個只有第一象限的極坐標，並將方差，相關係數進行捆綁，通過轉化為極坐標系坐標進行繪製。為了實現泰勒圖的繪製，我設計了兩個函數：

set_tayloraxes(fig, location=111) 和plot_taylor(axes, refsample, sample, args, *kwargs)

set_tayloraxes()函數用於建立一個泰勒圖的坐標系，這個自定義函數一般情況下不建議修改，每一個參數都是經過多次調試得到的，很可能牽一髮動全身。因此，將繪圖部分的獨立成為了plot_taylor函數()，這部分函數較為簡單，目的就是將需要繪圖的數據，轉換為極坐標系坐標，通過plot函數將散點打在泰勒圖上，這個函數模塊較為簡單，可以根據自己的輸入數據情況進行調整。

下面介紹下函數的具體用法：

輸入：

fig： 需要繪圖的figure

rect：圖的位置，如111為1行1列第一個，122為1行2列第2個

輸出：

polar_ax：泰勒坐標系

輸入：

axes : setup_axes返回的泰勒坐標系

refsample ：參照樣本

sample ：評估樣本

args, *kwargs ：plt.plot()函數的相關參數，設置點的顏色，形狀等等。

下面給出示例：

曼哈頓圖和QQ圖是兩個在全基因組關聯（GWAS）分析裡面最常出現的圖形，基本上已經是GWAS的標配，幾乎在每篇GWAS的文章都會見到，它們的作用和所要傳達出來的信息我也在上一篇關於GWAS的博文中做了些說明，在這裡我們就只集中在如何用Python和geneview將其有效地展現出來。

首先，準備一些數據來作為例子。

我這裡用來展現的數據是2011年丹麥人所做過的一個關於年輕人過度肥胖的GWAS研究——GOYA，數據也是從他們所發表的結果中獲得，總共有5,373個樣本，其中超重的個體（case）有2,633個，正常的個體（control）是2,740個，從樣本量上看還算可以。為了方便使用，我對其做了相關的處理，包括從PED和MAP文件到GEN文件的生成，並重複了一次case-control的關聯性分析，計算出了晶元上所研究的各個SNP位點與肥胖相關的顯著性程度（即p-value），最後又將結果數據抽取出來做成數據集——放在這裡供下載（15.6Mb，csv格式）。

【注】以上內容雖提及到了一些領域內術語和相關文件格式，但若不懂也請不必糾結，因為後續處理都是基於這個最終的數據集來完成的。

接著，需要將geneview軟體包加入到你的Python中，有多種不同的方式，但推薦直接使用pip，以下是安裝比較穩定的發布版，直接在終端命令行下(Linux or Mac)輸入：

pip install geneview

第三種辦法就是直接下載源碼，然後自行編譯，雖然不推薦這種做法（因為還有依賴包必須自行下載安裝，過程會比較麻煩低效），但對於某些不能連接外網的集群也只能如此，這三種方式都是可行的。

曼哈頓圖

將示例數據下載下來：

wget data/master/GOYA.csv

先簡單地查看一下數據的格式:

chrID,rsID,position,pvalue

1,rs3094315,742429,0.144586

1,rs3115860,743268,0.230022

1,rs12562034,758311,0.644366

1,rs12124819,766409,0.146269

1,rs4475691,836671,0.458197

1,rs28705211,890368,0.362731

1,rs13303118,908247,0.22912

1,rs9777703,918699,0.37948

1,rs3121567,933331,0.440824

一共是4列（逗號分隔），分別為：[1]染色體編號，[2]SNP rs 編號，[3] 位點在染色體上的位置，[4]顯著性差異程度（pvalue）。在本例曼哈頓圖中我們只需要使用第1,3和4列；而QQ圖則只需要第4列——pvalue。

下面先從繪製曼哈頓圖開始。我們先將需要的數據讀取到一個列表中，可以這樣做：

import csv

data = []

with open(“GOYA.csv”) as f:

f_csv = csv.reader(f)

headers = next(f_csv)

data = [[row[0], int(row[2]), float(row[3])] for row in f_csv]

現在GOYA.csv中的數據就都存放在data列表中了，由於Python在讀取文件中數據時，都是以string類型存放，因此對於第3和第4列的數據有必要事先把做點類型轉換。

接下來，調用geneview中的曼哈頓圖函數。

import matplotlib.pyplot as plt

from geneview.gwas import manhattanplot

ax = manhattanplot(data, xlabel=”Chromosome”, ylabel=”-Log10(P-value)”) # 這就是Manhattan plot的函數

plt.show()

只需這樣的一句代碼就能創建一個漂亮的曼哈頓圖，有必要再次指出的是，geneview是以matplotlib為基礎開發出來的，所創建的圖形對象實際上仍屬於matplotlib，geneview內部自定義了很多圖形風格，同時封裝了大量只屬於基因組數據的圖表類型，但圖形的輸出格式以及界面顯示都仍和matplotlib一樣，因此在這裡我們使用matplotlib.pyplot的show()函數(上例中：plt.show())將所繪製出來的曼哈頓圖顯示出來。如果要將圖形保存下來，則只需執行`plt.savefig(“man.png”)`，這樣就會在該目錄下生成一個名為『man.png』png格式的曼哈頓圖，若是要存為pdf格式，則只需將所要保存的文件名後綴改成『.pdf』（plt.savefig(“man.pdf”)）就可以了。下面這些格式：emf,

eps, pdf, png, jpg, ps, raw, rgba, svg,

svgz等都是支持的，至於最新的還有多少種，還請參照matplotlib文檔中說明。

此外，geneview中的每個畫圖函數都有著足夠的靈活性，我們也可以根據自己的需要做一些調整，比如：

xtick = [‘1’, ‘2’,’3′,’4′,’5′,’6′,’7′,’8′,’9′,’10’,’11’,’12’,’13’,’14’,’16’,’18’, ’20’,’22’]

manhattanplot(data,

xlabel=”Chromosome”, # 設置x軸名字

ylabel=”-Log10(P-value)”, # 設置y軸名字

xtick_label_set = set(xtick), # 限定橫坐標軸上的刻度顯示

s=40, # 設置圖中散點的大小

alpha=0.5, # 調整散點透明度

color=”#f28b1e,#9a0dea,#ea0dcc,#63b8ff”, # 設置新的顏色組合

)

實現新的顏色組合、限定x軸上的刻度顯示和散點大小的調節。甚至還可以將散點改為線：

manhattanplot(data，

xlabel=”Chromosome”, # 設置x軸名字

ylabel=”-Log10(P-value)”, # 設置y軸名字

xtick_label_set = set(xtick), # 限定橫坐標軸上的刻度顯示

alpha=0.5, # 調整散點透明度

color=”#f28b1e,#9a0dea,#ea0dcc,#63b8ff”, # 設置新的顏色組合

kind=”line”

)

其它方面的調整請查看geneview文檔中的相關說明。

Q-Q圖

qq圖只需用到上例中的pvalue那一列：

import csv

import matplotlib.pyplot as plt

from geneview.gwas import qqplot

pvalue=[]

with open(“GOYA.csv”) as f:

f_csv = csv.reader(f)

headers = next(f_csv)

pvalue = [float(row[3]) for row in f_csv]

ax = qqplot(pvalue, color=”#00bb33″, xlabel=”Expected p-value(-log10)”, ylabel=”Observed p-value(-log10)”) # Q-Q 圖

plt.show()

同樣，也可以根據自己的需要對改圖進行相關的調整。

以上，便是如何使用Python來製作Manhattan圖和QQ圖的方法，geneview的集成函數簡化了這樣的一個過程。

對於氣象繪圖來講，第一步是對數據的處理，通過各類公式，或者統計方法將原始數據處理為目標數據。

按照氣象統計課程的內容，我給出了一些常用到的統計方法的對應函數：

在計算氣候態，區域平均時均要使用到求均值函數，對應NCL中的dim_average函數，在python中通常使用np.mean()函數

numpy.mean(a, axis, dtype)

假設a為[time,lat,lon]的數據，那麼

需要特別注意的是，氣象數據中常有缺測，在NCL中，使用求均值函數會自動略過，而在python中，當任意一數與缺測(np.nan)計算的結果均為np.nan，比如求[1,2,3,4，np.nan]的平均值，結果為np.nan

因此，當數據存在缺測數據時，通常使用np.nanmean()函數，用法同上，此時[1,2,3,4，np.nan]的平均值為(1+2+3+4)/4 = 2.5

同樣的，求某數組最大最小值時也有np.nanmax(), np.nanmin()函數來補充np.max(), np.min()的不足。

其他很多np的計算函數也可以通過在前邊加『nan』來使用。

另外，

也可以直接將a中缺失值全部填充為0。

np.std(a, axis, dtype)

用法同np.mean()

在NCL中有直接求數據標準化的函數dim_standardize()

其實也就是一行的事，根據需要指定維度即可。

皮爾遜相關係數：

相關可以說是氣象科研中最常用的方法之一了，numpy函數中的np.corrcoef(x, y)就可以實現相關計算。但是在這裡我推薦scipy.stats中的函數來計算相關係數：

這個函數缺點和有點都很明顯，優點是可以直接返回相關係數R及其P值，這避免了我們進一步計算置信度。而缺點則是該函數只支持兩個一維數組的計算，也就是說當我們需要計算一個場和一個序列的相關時，我們需要循環來實現。

其中a[time,lat,lon]，b[time]

(NCL中為regcoef()函數)

同樣推薦Scipy庫中的stats.linregress(x,y)函數：

slop: 回歸斜率

intercept：回歸截距

r_value： 相關係數

p_value： P值

std_err： 估計標準誤差

直接可以輸出P值，同樣省去了做置信度檢驗的過程，遺憾的是仍需同相關係數一樣循環計算。