python将等值线数据转二维标量场(曲面重采样)

将等值线转二维标量场contour2xyz

二维标量场的表示方式有很多，这里介绍将其等值线数据转为二维平面等间隔采样的数据。

曲面的两种采样方式

已知等值线的经纬度数据（z等间隔），想将它重采样为经纬度等间隔的二维曲面（xy等间隔）。
数据文件格式为：

> contour 25 
-137 45 25
-138 45 25
-139 42 25
> contour 50
-134 44 50
-135 45 50
...
> contour 75
...

提取曲面轮廓

使用scipy.spatial.ConvexHull提取曲面经纬度范围。

from scipy.spatial import ConvexHull
import matplotlib.pyplot as plt
import numpy as np 

#读入等值线的xyz数据
ctr1 = np.loadtxt("contour.dat",comments=">") 

#提取外部范围
hull1 = ConvexHull( points=ctr1[:,[0,1]] ) 
#范围的经纬度曲线 
bndy1 = np.vstack((ctr1[hull1.vertices,0],ctr1[hull1.vertices,1])).T 

#画外部范围
#方法1
plt.plot( x=ctr1[hull1.vertices,0], y=ctr1[hull1.vertices,1] ) 
#方法2
for i in hull1.simplices:
    fig.plot(x=ctr1[i,0],y=ctr1[i,1])
plt.show()

制作轮廓掩膜

这一步骤主要是为了防止对曲面以外的数据进行采样，所以制作一个mask

#确定掩膜mask的经纬度间隔
out = []
lon = np.arange(-180,180+0.1,0.1)
lat = np.arange(-90, 90+0.1, 0.1)

#判断点(xx,yy)是否在轮廓线bndy1的范围内（射线法）
#函数isin_multipolygon拷贝自https://www.cnblogs.com/shld/p/9758303.html
x,y = np.meshgrid(lon,lat)
for xx,yy in np.vstack((x.flatten(), y.flatten())).T:
    zz = 1*isin_multipolygon( [xx,yy], bndy1.tolist(), contain_boundary=True) 
    out.append([xx,yy,zz])

#生成掩膜mask
msk = np.array(out)

曲面重采样

使用scipy.interpolate.griddata对等值面重采样

from scipy.interpolate import griddata 

#重采样的经纬度范围
lon = np.arange( -180, 180+0.1, 0.1 )
lat = np.arange( -90,   90+0.1, 0.1 )

#设置重采样等间距网格
x,y = np.meshgrid( lon, lat )
loc = np.vstack(( x.flatten(), y.flatten() )).T

#重采样并保存文件
h   = griddata( points=ctr1[:,[0,1]], values=ctr1[:,2], xi=loc )
out = np.vstack(( loc.T, h*msk[:,2] )).T
np.savetxt("contour.xyz",out,fmt="%6.1f %4.1f %6.2f",header="lon. lat. value")

输出结果

-180.0 -90.0 NaN
...
-134.4 34.4 23.57
-134.5 34.4 22.46
...