python寫入二維netCDF文件（txt轉nc文件）

netCDF文件后綴名是.nc，是一種二進制的網格文件，無法直接打開，下面給出基於python的netCDF文件生成的腳本。

（該程序目前已發布更新版本，鏈接是https://www.cnblogs.com/liangxuran/p/13659283.html）

輸入數據說明：

數據為31°N-40°N，130°E-142°E區域表層的地震波速，網格分布大小不均勻。數據存儲在txt文件中，分為三列，無表頭。

部分數據記錄見下圖，第一列是經度，第二列是維度，第三列是波速擾動值%

 

過程與代碼 

為了將其生成均勻網格的二維nc文件，我們需要1.文件讀取，2.數據插值，3.寫入nc文件。

'''
program: ncdata.py
goal:    interpolate txt data and generate grd file
input:   vp1_02.txt
output:  vp1_03.nc
author:  Liang Xuran
e-mail:  [e-mail]
time:    2020.08.23
'''
import netCDF4 as nc
import pandas as pd
import numpy as np
import copy

file_txt = "vp1_02.txt"#txt文件位置
file_nc  = "vp1_03.nc" #nc文件輸出位置
npa = 37  #緯向格點數(uninpolate)
nra = 45  #經向格點數(uninpolate)
npx = 180 #緯向格點數(inpolate)
nrx = 240 #經向格點數(inpolate)

#使用Dataset創建nc文件
nc_pro = nc.Dataset("vp1_03.nc", 'w', format = 'NETCDF4')

#命名nc文件變量，None為自由維度,name='pna',size=npa
nc_pro.createDimension('pna',npx+1)
nc_pro.createDimension('rna',nrx+1)

#設置nc文件中變量類型，這里使用的是float64
pna = nc_pro.createVariable('pna',np.float64,('pna',))
rna = nc_pro.createVariable('rna',np.float64,('rna',))
dv  = nc_pro.createVariable('dv',np.float64,('pna','rna'))

#定義變量單位
pna.units = 'degree'
rna.units = 'degree'
dv.units  = '%'

#使用panda中的read_csv讀取txt文件
#type(nc_data)=<class 'pandas.core.frame.DataFrame'>
nc_data=pd.read_csv("vp1_02.txt",delim_whitespace=True,names='rpv')

#插值操作:精度0.05degree,范圍rna=130:142;pna=31:40
pnx = nc_data['p'].values[0:1665:nra]
rnx = nc_data['r'].values[0:nra]
dvx = nc_data['v'].values.reshape(npa,nra)
rx  = np.linspace(130,142,241)
px  = np.linspace(31,40,181)
vx  = np.full([181,241],np.nan)
for i in range(npa):
    for j in range(nra):
        m = int(round((pnx[i]-pnx[0])/0.05))
        n = int(round((rnx[j]-rnx[0])/0.05))
        vx[m,n]=dvx[i,j]
#縱向插值,從vx[181,241]獲得中間量dv_inpolate1[181,241]
nc_v_uninpo  = pd.DataFrame(vx) #transfer np.array into pandas.DataFrame
nc_v_inpo1   = nc_v_uninpo.interpolate() #pandas.DataFrame interpolate
dv_inpolate1 = nc_v_inpo1.values #tansfer pandas.DataFrame into np.array
#橫向插值,從dv_inpolate1[181,241]獲得中間量dv_inpolate2[181,241]
nc_v_inpo    = pd.DataFrame(dv_inpolate1.T)
nc_v_inpo2   = nc_v_inpo.interpolate()
dv_inpolate2 = nc_v_inpo2.values.T

#將數值輸入到nc文件中,input must be np.array
nc_pro.variables['pna'][:]  = px
nc_pro.variables['rna'][:]  = rx
nc_pro.variables['dv'][:,:] = dv_inpolate2[:,:]
#為避免內存的占用，寫完一個關閉一個
nc_pro.close()
    

 

踩坑指南：

1.文件路徑最好寫相對路徑，第15-16行寫絕對路徑/mnt/f/research/vp1_02.txt會出現報錯。

2.插值前后格點數量有所擴充，由於經緯度的變量檢索是從0開始的，因此，26-27行創建變量的維度需要加1，也可以定位None，代表自由維度。

3.元組只有一個元素時，需要在括號末尾加“，”，例如30-31行。

4.用pandas.read_csv讀取文件時，第41行，默認數據之間時用“，”逗號連接的，需要修改默認值delim_whitespace代表用空格連接。names代表每一列的元素名稱。

5.二維數據pandas.DataFrame進行賦值/輸出內容時，需要使用其values屬性，例如44-46行。在插值過程中，也需要將numpy.array矩陣轉為pandas.DataFrame網格后進行，例如55-62行。

6.nc文件的變量被賦值的時候，只接受numpy.array數據，例如65-67行；二維的numpy.array數據索引是一個中括號而不是兩個中括號[0:181,0:241]，例如54行和67行。

7.python的浮點型數據和整形數據加減運算的時候要小心，例如52-53行計算mn的時候，出現了以下現象：

8.關於變量的索引和插值前后數據的對比檢驗。可以參考下圖：

 

檢驗未插值部分的數據，選擇34.25°N-34.30°N，134.80°E-134.85°E區域： dvx[10:12,12:14]  vx[65:67,96:98]

檢驗已插值部分的數據，選擇33.2°N-33.4°N，136.6°E-136.8°E的區域： dvx[2:4,41:43]  vx[4:9,132:137]