python 创建netcdf_python写入二维netCDF文件进阶版程序（txt转nc文件）

之前发布了用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.py3 goal: generate grd file4 input: vp1_nit1.txt5 output: vp1_CD.nc6 author: Liang Xuran7 e-mail: xuranliang@gig.ac.cn8 time: 2020.09.109 '''

10 importnetCDF4 as nc11 importpandas as pd12 importnumpy as np13 importmath14

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 defreadfile(file_in,npa,nra,nha):33 #使用panda中的read_csv读取txt文件

34 #type(pd_data)=

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).T41 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 defpointloc(pnx,rnx,hnx,x):48 #判断点x在第几个格点之间

49 #x的左侧/右侧/下侧/上侧第几个网格线

50 left = -1 + sum( i<= x[0] for i inrnx )51 right = rnx.size - sum( i>= x[0] for i inrnx )52 down = -1 + sum( i<= x[1] for i inpnx )53 up = pnx.size - sum( i>= x[1] for i inpnx )54 top = -1 + sum( i<= x[2] for i inhnx )55 bottom= hnx.size - sum( i>= x[2] for i inhnx )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 defpointwv(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 *lenh84 #经线分成上下两段[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 /volumn104 wv[1] = lenp1 * lenr1 * lenh0 /volumn105 wv[2] = lenp0 * lenr0 * lenh0 /volumn106 wv[3] = lenp1 * lenr0 * lenh0 /volumn107 wv[4] = lenp0 * lenr1 * lenh1 /volumn108 wv[5] = lenp1 * lenr1 * lenh1 /volumn109 wv[6] = lenp0 * lenr0 * lenh1 /volumn110 wv[7] = lenp1 * lenr0 * lenh1 /volumn111 returnnp.array(wv)112

113 defpointv(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 returnnp.array(v)126

127 defwritefile(filename,din,hin,vin):128 nc_pro = nc.Dataset(filename, 'w', format='NETCDF4')129 lenx =din.size130 leny =hin.size131 #命名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'][:] =hin144 nc_pro.variables['d'][:] =din145 nc_pro.variables['v'][:,:] =vin146 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 inrange(nody):162 for i inrange(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网格文件函数保留了上一版本的内容，相关参考上一篇博客。