之前發布了用python將txt轉netCDF格式的文章（鏈接是https://www.cnblogs.com/liangxuran/p/13551134.html），網上被轉的有很多。但是那個是最初的版本，有以下幾個缺點：
1.插值過程是硬傷，pandas自帶的插值函數是一維插值，因此在進行三維插值時需要改變三次方向，最終插值的結果並不理想。
2.只能生成某深度/經度/緯度剖面的二維的nc文件，無法生成已知起始點和終點經緯度的測線剖面。
經過我的多次改良，發布了v5.0版本，ncdata_v5.0.py具有以下優良性質：
1.使用三維線性插值，空間某點的函數值由周圍的8個格點（grid）決定，距離越近，權重越大。
2.快速，插值函數由幾個函數確定，沒有用到pandas自帶的插值函數，運行速度提升一倍。
3.支持任何方向的側線，只需要輸入起點和終點的經緯度，就可以生成該側線的剖面函數結構（插值后的結果）
4.自定義分辨率，自定義netCDF文件的橫縱坐標網格數量。
另外，關於網格設置，通常網格分為兩種，一種grid型，一種cell型。grid型網格是將空間中特定點的函數值計算出來，其余位置的函數值由周圍網格插值出來。cell型網格則是將空間分成若干個晶胞，每個晶胞內部的函數值是相同的，即晶胞內部均勻，晶胞之間有差異。本程序使用的是grid型網格。
關於程序的輸入和輸出文件，上一篇博客有詳細的說明，如果閱讀過程有不理解的，可以先看上一篇文章（鏈接是https://www.cnblogs.com/liangxuran/p/13551134.html）
1 ''' 2 program: ncdata_v5.py 3 goal: generate grd file 4 input: vp1_nit1.txt 5 output: vp1_CD.nc 6 author: Liang Xuran 7 e-mail: xuranliang@gig.ac.cn 8 time: 2020.09.10 9 ''' 10 import netCDF4 as nc 11 import pandas as pd 12 import numpy as np 13 import math 14 15 #------------運行前修改以下內容------------ 16 file_txt= "vp1_nit1.txt" #輸入txt文件 17 file_nc = "vp1_CD.nc" #輸出nc文件名 18 x_a = [134.80 , 34.80] #起點坐標經緯度 19 x_b = [135.30 , 34.40] #終點坐標經緯度 20 npa = 37 #緯向格點數(uninpolate) 21 nra = 45 #經向格點數(uninpolate) 22 nha = 10 #徑向格點數(uninpolate) 23 ranx = 63.82 #新剖面x軸范圍 24 rany = 25 #新剖面y軸范圍 25 nodx = 40 #側線插值點數 26 nody = 26 #深度插值點數 27 #注意ranx即剖面測線長度，可用球面距離公式/GMT軟件求出 28 #R0·arccos[cospA*cospB*cos(rA-rB)+sinpA*sinpB] 29 #echo '$B' | gmt project -C$A -E$B -Fpz -Q | cat 30 #------------修改完畢后開始運行------------ 31 32 def readfile(file_in,npa,nra,nha): 33 #使用panda中的read_csv讀取txt文件 34 #type(pd_data)=<class 'pandas.core.frame.DataFrame'> 35 pd_data=pd.read_csv(file_in,delim_whitespace=True,names='prhv') 36 #粗略3d網格生成,平面數據轉變為3d非均勻立體網格數據 37 pnx = pd_data['p'].values[0:nra*npa:nra] 38 rnx = pd_data['r'].values[0:nra] 39 hnx = pd_data['h'].values[0:nra*npa*nha:nra*npa] 40 dvx1= pd_data['v'].values.reshape(nha,npa*nra).T 41 dvx = dvx1[:,0].reshape(npa,nra) 42 for k in range(1,nha): 43 dvx2 = dvx1[:,k].reshape(npa,nra) 44 dvx = np.dstack((dvx,dvx2)) 45 return [pnx,rnx,hnx,dvx] 46 47 def pointloc(pnx,rnx,hnx,x): 48 #判斷點x在第幾個格點之間 49 #x的左側/右側/下側/上側第幾個網格線 50 left = -1 + sum( i<= x[0] for i in rnx ) 51 right = rnx.size - sum( i>= x[0] for i in rnx ) 52 down = -1 + sum( i<= x[1] for i in pnx ) 53 up = pnx.size - sum( i>= x[1] for i in pnx ) 54 top = -1 + sum( i<= x[2] for i in hnx ) 55 bottom= hnx.size - sum( i>= x[2] for i in hnx ) 56 if left == rnx.size-1: left = left - 1 57 if right == 0 : right = right + 1 58 if down == pnx.size-1: down = down - 1 59 if up == 0 : up = up + 1 60 if top == hnx.size-1: top = top - 1 61 if bottom== 0 : bottom= bottom+ 1 62 if left == right : right = right + 1 63 if down == up : up = up + 1 64 if top == bottom : bottom= bottom+ 1 65 leupbo = [ rnx[left] , pnx[up] , hnx[bottom] ] 66 ledobo = [ rnx[left] , pnx[down] , hnx[bottom] ] 67 riupbo = [ rnx[right] , pnx[up] , hnx[bottom] ] 68 ridobo = [ rnx[right] , pnx[down] , hnx[bottom] ] 69 leupto = [ rnx[left] , pnx[up] , hnx[top] ] 70 ledoto = [ rnx[left] , pnx[down] , hnx[top] ] 71 riupto = [ rnx[right] , pnx[up] , hnx[top] ] 72 ridoto = [ rnx[right] , pnx[down] , hnx[top] ] 73 pos8=[leupbo,ledobo,riupbo,ridobo,leupto,ledoto,riupto,ridoto] 74 loc6 = [ left, right, down, up, bottom, top ] 75 return [pos8,loc6] 76 77 def pointwv(pos8,x): 78 #計算網格點的權重 79 #大網格長寬(右上底點2-左下頂點5)(0-x;1-y;2-z) 80 lenr = pos8[2][0]-pos8[5][0] 81 lenp = pos8[2][1]-pos8[5][1] 82 lenh = pos8[2][2]-pos8[5][2] 83 volumn= lenr * lenp * lenh 84 #經線分成上下兩段[1]-r方向 85 lenr1 = pos8[2][0]-x[0] #下段長，右上底點2-x點 86 lenr0 = x[0]-pos8[5][0] #上段長，x點-左下頂點5 87 #緯線分成左右兩段[0]-p方向 88 lenp1 = pos8[2][1]-x[1] #左段長，右上底點2-x點 89 lenp0 = x[1]-pos8[5][1] #右段長，x點-左下頂點5 90 #深度線分成深淺兩段[2]-h方向 91 lenh1 = pos8[2][2]-x[2] #淺段長，右上底點2-x點 92 lenh0 = x[2]-pos8[5][2] #深段長，x點-左下頂點5 93 #權重分配,等效於頂點對角小立方體的體積百分比 94 # p leny0+leny1 95 # leupto 4|_______6 riupto 96 # /| | 97 # leupbo 0/ | | 98 # | |_______|___r lenx0+lenx1 99 # | / 5 / 7 ridoto 100 # ledobo 1 |/_______/ 3 ridobo 101 # /h lenz0+lenz1 102 wv = [0,0,0,0,0,0,0,0] 103 wv[0] = lenp0 * lenr1 * lenh0 / volumn 104 wv[1] = lenp1 * lenr1 * lenh0 / volumn 105 wv[2] = lenp0 * lenr0 * lenh0 / volumn 106 wv[3] = lenp1 * lenr0 * lenh0 / volumn 107 wv[4] = lenp0 * lenr1 * lenh1 / volumn 108 wv[5] = lenp1 * lenr1 * lenh1 / volumn 109 wv[6] = lenp0 * lenr0 * lenh1 / volumn 110 wv[7] = lenp1 * lenr0 * lenh1 / volumn 111 return np.array(wv) 112 113 def pointv(loc,dvx): 114 #計算中心網格點的速度 115 #loc是一維np數組0左,1右,2下,3上,4底,5頂 116 v = [0,0,0,0,0,0,0,0] 117 v[0] = dvx [ loc[0] , loc[3] , loc[4] ] 118 v[1] = dvx [ loc[0] , loc[2] , loc[4] ] 119 v[2] = dvx [ loc[1] , loc[3] , loc[4] ] 120 v[3] = dvx [ loc[1] , loc[2] , loc[4] ] 121 v[4] = dvx [ loc[0] , loc[3] , loc[5] ] 122 v[5] = dvx [ loc[0] , loc[2] , loc[5] ] 123 v[6] = dvx [ loc[1] , loc[3] , loc[5] ] 124 v[7] = dvx [ loc[1] , loc[2] , loc[5] ] 125 return np.array(v) 126 127 def writefile(filename,din,hin,vin): 128 nc_pro = nc.Dataset(filename, 'w', format='NETCDF4') 129 lenx = din.size 130 leny = hin.size 131 #命名nc文件變量，None為自由維度,name='hna',size=nha 132 nc_pro.createDimension('h',leny) 133 nc_pro.createDimension('d',lenx) 134 #設置nc文件中變量類型，這里使用的是float64 135 h = nc_pro.createVariable('h', np.float32, ('h',) ) 136 d = nc_pro.createVariable('d', np.float32, ('d',) ) 137 v = nc_pro.createVariable('v', np.float32, ('h','d')) 138 #定義變量單位 139 h.units = 'km' 140 d.units = 'km' 141 v.units = 'percentage' 142 #將數值輸入到nc文件中,input must be np.array 143 nc_pro.variables['h'][:] = hin 144 nc_pro.variables['d'][:] = din 145 nc_pro.variables['v'][:,:] = vin 146 print("nc file messeges:",nc_pro) 147 #為避免內存的占用，寫完一個關閉一個 148 nc_pro.close() 149 150 #主程序開始 151 [pna,rna,hna,dva]= readfile (file_txt,npa,nra,nha) 152 #xab[i,0:1]代表測線AB上的第i個插值點的經緯 153 xab = np.zeros ( (nodx , 2) ) 154 xab[:,0] = np.linspace ( x_a[0], x_b[0], nodx ) 155 xab[:,1] = np.linspace ( x_a[1], x_b[1], nodx ) 156 #插值平面的xy坐標值，用於給nc文件賦值，單位km 157 corx = np.linspace ( 0, ranx, nodx ) 158 cory = np.linspace ( 0, rany, nody ) 159 #插值平面的網格數值 160 vx = np.zeros ( (nody,nodx) ) 161 for h in range(nody): 162 for i in range(nodx): 163 #第h層第i個插值點的經緯深坐標是 164 #經度xab[i,0]緯度xab[i,1]深度cory[h] 165 #中心點的坐標 166 cor = np.array([ xab[i,0], xab[i,1], cory[h] ]) 167 #周圍8個頂點的坐標 168 [pos,loc] = pointloc(pna,rna,hna,cor) 169 #周圍8個頂點的權重 170 wv8 = pointwv (pos,cor) 171 #周圍8個頂點的速度 172 v8 = pointv (loc,dva) 173 #核心點的速度 174 vx[h,i] = np.dot(v8.T , wv8) 175 #將插值平面的網格數值寫入nc文件 176 writefile( file_nc, corx, cory, vx)
寫代碼過程有點辛苦，debug一晚上終於成功了，踩坑指南有很多，一個一個函數的說：
（1）第41-44行讀取函數時，txt中的函數數組需要從一維變為三維，在reshape的時候，先按深度分為nha層，然后每層分別reshape和stack，最終生成三維函數數組dvx。
（2）第50-64行判斷某點位於哪些格點之間時，首先需要考慮三種情況：（1）該點位於網格下界上，屬於該網格，（2）該點位於網格上界上，屬於下一網格，（3）該點位於網格上下界之間，屬於該網格。然后還需要考慮兩種情況：（4）該點位於整個區域下界，屬於第一個網格。（5）該點位於整個區域上界，屬於最后一個網格。
（3）第65-74行判斷某點位置時，該程序輸出的是該網格周圍是8個頂點和6條棱。其中網格簡圖在93-101行畫出來了。
（4）第80-92行在空間線性插值的時候，需要知道該點在該網格中的具體位置，也就是過該點的三個方向的面把立方體分成8個小立方體，即該點把網格的長寬高都分成兩段，這6段長度只需要知道三個點就可以算出來：插值點的坐標（經緯深度），右上底點坐標（即93-101行簡圖中的2點），坐下頂點（即93-101簡圖中的5點）。
（5）第102-110行空間插值的權重分配。某網格某頂點的權重=插值點與對角點形成的小立方體體積/該網格體積
（6）寫入netCDF網格文件函數保留了上一版本的內容，相關參考上一篇博客。