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

 

posted @ 2020-09-24 10:10  Philbert  阅读(1644)  评论(1)  编辑  收藏  举报