之前发布了用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网格文件函数保留了上一版本的内容，相关参考上一篇博客。