海洋数据处理与插值
0.数据处理和形式转换
0.1.数据处理
0.1.1.HYCOM数据处理
1.插值
"""
1.数据介绍
HYCOM数据是低分辨率数据，下面用到的分辨率是0.08×0.04(lon, lat)，深度已经将其缩减到0-800m的26个不等深度层：
depth_l = [0, 10, 15, 20, 25, 30, 35, 40, 45, 50, 60, 70, 80, 90, 1000, 125, 150, 200, 250, 300, 350, 400, 500, 600, 700, 800]
经过处理后的数据形式有两种，方便不同的擦插值方法和画图方法：
1）：
lat lon 0 10 ... 500 600 700 800
0 20.52 118.96 28.449001 28.472000 ... 8.603999 7.384000 6.513000 5.709999
1 20.52 119.04 28.520000 28.471001 ... 8.556000 7.419999 6.514999 5.717999
2 20.52 119.12 28.472000 28.483002 ... 8.544000 7.384000 6.542999 5.732999
3 20.52 119.20 28.479000 28.488001 ... 8.493999 7.398999 6.565999 5.744999
4 20.52 119.28 28.490002 28.514999 ... 8.474999 7.412999 6.563999 5.728999
5 20.52 119.36 28.504002 28.530001 ... 8.377000 7.457999 6.549000 5.742999
2):就是把第一种的深度数据堆叠形成新的一列，这里和我们glider数据形式相同，可以方便混合在一起进行插值
lat lon temp depth
0 20.52 118.96 28.449001 0
1 20.52 119.04 28.520000 0
2 20.52 119.12 28.472000 0
3 20.52 119.20 28.479000 0
4 20.52 119.28 28.490002 0
"""
1.1.三维插值
1.1.1.用pykring进行三维克里金插值
from pykrige.ok3d import OrdinaryKriging3D
from pykrige.ok import OrdinaryKriging
# 准备插值后的网格数据，这里的hycom数据是第二种形式的数据
grid_lon = np.linspace(hycom['lon'].values[0], hycom['lon'].values[-1], 100)
grid_lat = np.linspace(hycom['lat'].values[0], hycom['lat'].values[-1], 100)
grid_depth = np.linspace(0, 800, 8) # 20m一个深度
# 插值之前的轴坐标，和对应数据
lons = hycom['lon'].values
lats = hycom['lat'].values
depth = hycom['depth'].values
# 普通克里金3维插值类
OK = OrdinaryKriging3D(lons, lats, depth, hycom['temp'].values, variogram_model='gaussian')
z1, ss1 = OK.execute('grid', grid_lon, grid_lat, grid_depth)
print(z1.shape) # z1就是插值后的网格数据,形状是（8， 100， 100）
注意上述程序因为插值后的网格点太多，我没跑出来结果，应该是时间不够，因为我把插值后的网格改少之后是可以跑通的。比如lat和lon上是10， 深度上是3层
总之，这个不知道啥原因，三维插值很慢很慢，数据大了根本用不了
1.1.2.用scipy的griddata方法进行三维插值，其中关键字为nearst和linear时才能进行三维插值，cubic不能
from scipy.interpolate import griddata
# 准备插值后的网格数据
grid_lon = np.linspace(hycom['lon'].values[0], hycom['lon'].values[-1], 100)
grid_lat = np.linspace(hycom['lat'].values[0], hycom['lat'].values[-1], 100)
grid_depth = np.linspace(0, 800, 8) # 20m一个深度
x, y, z = np.meshgrid(grid_lon, grid_lat, grid_depth)
data_inter = griddata(hycom[['lon', 'lat', 'depth']].values, hycom['temp'].values, (x, y, z), method='linear')
print(data_inter.shape) # (100, 100, 8)
这个速度很快，不过效果略差，感觉不如cubic和克里格插值法，不过目前也只能采用这个了
1.2.二维插值
1.2.1.二维克里金插值
# 准备插值后的网格数据
grid_lon = np.linspace(hycom['lon'].values[0], hycom['lon'].values[-1], 100)
grid_lat = np.linspace(hycom['lat'].values[0], hycom['lat'].values[-1], 100)
# 插值之前的轴坐标，和对应数据，这里的hycom_是第一种形式的数据
lons = hycom_['lon'].values
lats = hycom_['lat'].values
temp = hycom_[0].values # 0m处的温度值
# 二维普通克里金
OK = OrdinaryKriging(lons, lats, temp, variogram_model='linear')
z1, ss1 = OK.execute('grid', grid_lon, grid_lat)
print(z1.shape) # （100，100）
插值速度比cubic慢一些，不过也挺快的了，效果也很好，不过和cubic（三次样条插值-二维）肉眼看差距不大
1.2.2.scipy的griddata的cubic（三次样条插值）
from scipy.interpolate import griddata
lat_lon_points = np.array(hycom_[['lat', 'lon']], dtype=np.float32) # 查之前的数据点，注意这里不是网格点，是一个一个的（lat,lon）点
grid_x, grid_y = np.mgrid[20.53:21.52:100j, 119.01:120:100j] # 插值后的网格点
values = data[0].values # 0m处的温度值
values_inter = griddata(lat_lon_points, values, (grid_x, grid_y), method='cubic')
print(values_inter.shape) # (100, 100)
效果看着和克里格没啥区别，速度也是
1.2.3.scipy的griddata的linear二维插值
from scipy.interpolate import griddata
lat_lon_points = np.array(hycom_[['lat', 'lon']], dtype=np.float32) # 查之前的数据点，注意这里不是网格点，是一个一个的（lat,lon）点
grid_x, grid_y = np.mgrid[20.53:21.52:100j, 119.01:120:100j] # 插值后的网格点
values = data[0].values # 0m处的温度值
values_inter = griddata(lat_lon_points, values, (grid_x, grid_y), method='linear')
print(values_inter.shape) # (100, 100)
效果看起来也可以，不过感觉比cubic和克里金插值法稍微差一点
1.2.4.scipy的Rbf进行二维插值
from scipy.interpolate import Rbf
grid_lon = np.linspace(hycom['lon'].values[0], hycom['lon'].values[-1], 100)
grid_lat = np.linspace(hycom['lat'].values[0], hycom['lat'].values[-1], 100)
rbf = Rbf(hycom_['lon'], hycom_['lat'], hycom_[0])
XI, YI = np.meshgrid(grid_lon, grid_lat)
ZI = rbf(XI, YI)
print(ZI.shape) # （100， 100）
效果肉眼上介于cubic和线性插值之间吧
1.3.总结
"""
目前三维插值选griddata，方法选择'linear'。pykrige三维插值太慢，rbf三维插值说我的矩阵异常，结果是错的，插不出来。
二维就选择griddata的cubic或者克里金吧！
"""
2.画图
2.1.二维图
2.1.1.二维温度云图，plt.imshow()
def plot_hot(data, depth=0, save=False):
"""画不同深度的热力图，默认只画0m深度的，也可以传入一个包含四个元素的深度列表，要保存图片的话，save传入文件名"""
if isinstance(depth, list): # 如果要画四个深度层在一个图上，也就是画四个子图
for i, dt in enumerate(data[depth].values.T):
plt.subplot(2, 2, i + 1)
im = plt.imshow(dt.reshape(-1, 100), cmap=plt.cm.hot)
plt.xticks([0, 50, 99], ['119°E', '119.5°E', '120°E'], fontsize=7)
# 其实是21.52-22.52
plt.yticks([0, 25, 50, 75, 99], ['21.5°N', '21.75°N', '22.0°N', '22.25°N', '22.5°N'], fontsize=7)
plt.title(f'{depth[i]}m')
cbar = plt.colorbar(im)
cbar.ax.set_title('T(°C)', loc='left', fontsize=7) # 设置colorbar的标题
cbar.ax.tick_params(labelsize=7)
plt.subplots_adjust(wspace=0.3, hspace=0.3)
if save:
plt.savefig(f'./{save}_list.png', dpi=800, bbox_inches='tight')
plt.show()
else:
plt.figure()
im = plt.imshow(data[depth].values.reshape(-1, 100), cmap='hot')
plt.xticks([0, 50, 99], ['119°E', '119.5°E', '120°E'], fontsize=10)
# 其实是21.52-22.52
plt.yticks([0, 25, 50, 75, 99], ['21.5°N', '21.75°N', '22.0°N', '22.25°N', '22.5°N'], fontsize=10)
plt.title(f'{depth}m')
cbar = plt.colorbar(im)
cbar.ax.set_title('T(°C)', loc='left', fontsize=10) # 设置colorbar的标题
if save:
plt.savefig(f'./{save}.png', dpi=800, bbox_inches='tight')
plt.show()
2.1.2.二维等温线图，plt.contour()
def plot_contour(data, depth=0, save=False, num_col=100):
"""
# 根据传入的数据data：dataframe格式，形状是（10000， 28），后面要改进，可以处理不同形状的数据
画不同深度的温度等值线图，这里需要传入一个深度列表
"""
if isinstance(depth, list):
Y = data['lat'].values.reshape(-1, num_col)
X = data['lon'].values.reshape(-1, num_col)
fig, ax = plt.subplots(2, 2)
ax = ax.flatten()
for i, d in enumerate(depth):
Z = data[d].values.reshape(-1, num_col)
CS = ax[i].contour(X, Y, Z)
ax[i].clabel(CS, inline=True, fontsize=7)
ax[i].set_title(f'{d}m')
ax[i].get_xticklabels()
ax[i].set_xticks([119.2, 119.5, 119.8])
ax[i].set_xticklabels(['119.2°E', '119.5°E', '119.8°E'], fontsize=7)
ax[i].set_yticks([20.7, 21.0, 21.3])
ax[i].set_yticklabels(['20.7°N', '119.5°N', '119.8°N'], fontsize=7)
plt.subplots_adjust(wspace=0.3, hspace=0.3)
if save:
plt.savefig(f'./{save}_list.png', dpi=800, bbox_inches='tight')
plt.show()
else:
Y = data['lat'].values.reshape(-1, num_col)
X = data['lon'].values.reshape(-1, num_col)
plt.figure()
Z = data[depth].values.reshape(-1, num_col)
CS = plt.contour(X, Y, Z)
plt.clabel(CS, inline=True)
plt.title(f'{depth}m')
plt.xticks([119.2, 119.5, 119.8], ['119.2°E', '119.5°E', '119.8°E'])
plt.yticks([20.7, 21.0, 21.3], ['20.7°N', '119.5°N', '119.8°N'])
if save:
plt.savefig(f'./{save}.png', dpi=800, bbox_inches='tight')
plt.show()
2.1.3.二维散点图和折线图，plt.scatter(), 主要是设置grid来显示我们的数据在二维上的分布情况
def plot_scatter(data):
"""
:param data: 滑翔机数据，第二种形式的dataframe
:return: 展示滑翔机平面轨迹图
"""
point_list = extract_node(data)
lon = [x[0] for x in point_list]
lat = [x[1] for x in point_list]
plt.figure()
plt.scatter(lon[-20:], lat[-20:], s=12, c='red')
plt.plot(lon[-20:], lat[-20:], 'k--', alpha=0.5)
plt.show()
# 画数据融合示意图--散点图
fig = plt.figure()
ax = fig.add_subplot()
ax.scatter(sample_point_100['lon'], sample_point_100['lat'], label='采样点', s=3, c='red')
ax.scatter(hycom.data['lon'], hycom.data['lat'], label='预测点', s=4, c='blue')
maloc = plt.MultipleLocator(0.16)
miloc = plt.MultipleLocator(0.08)
ax.xaxis.set_major_locator(maloc)
ax.xaxis.set_minor_locator(miloc)
maloc_y = plt.MultipleLocator(0.04)
ax.yaxis.set_minor_locator(maloc_y)
ax.set_xlim(118.96, 120)
ax.set_ylim(20.52, 21.52)
# ax.set_xticks(np.arange(118.96, 120.01, 0.08), minor=True)
# ax.set_yticks(np.arange(20.52, 21.521, 0.04), minor=True)
ax.set_xlabel('Longitude/°E')
ax.set_ylabel('Latitude/°N')
ax.grid(which='both', axis='both', alpha=0.4)
plt.legend(loc='upper right')
# plt.savefig('./数据融合_100个采样点与364个预测点的分布示意图.png', dpi=800, bbox_inches='tight')
plt.show()
2.2.三维图
2.2.1.三维散点图和折线图，plt.scatter3D, plt.plot3D(), 主要显示glider每个剖面运动情况
def plot_scatter3D(data, profile):
plt.figure()
ax = plt.axes(projection='3d')
data = data[data['profile']==profile] # profile是剖面号
x, y, z = data['lon'], data['lat'],data['depth']
ax.scatter3D(x, y, z, alpha=0.5)
ax.plot3D(x, y, z, c='red')
ax.invert_zaxis() # 反转z轴，让z轴数值从大到小
ax.set(
xlabel = 'Longitude/°E',
ylabel = 'Latitude/°N',
zlabel = 'Depth/m',
title = f'剖面{profile}轨迹图',
xticks = [data['lon'].min(), data['lon'].max()],
yticks = [data['lat'].min(), data['lat'].max()],
xticklabels = [round(data['lon'].min(), 3), round(data['lon'].max(), 3)],
yticklabels = [round(data['lat'].min(), 3), round(data['lat'].max(), 3)]
# xlim = (120.04, 120.111),
# ylim = (20.695, 21.707),
)
# ax.view_init(180, 90)
plt.savefig('./1.png', dpi=800, bbox_inches='tight')
plt.show()
2.2.2.三维等温线图，显示表面0m处海水大致情况，垂直剖面海水变化情况，contourf（）
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
x = np.arange(118.96, 120.01, 0.08)
y = np.arange(20.52, 21.521, 0.04)
z = depth_l
X, Y, = np.meshgrid(x, y) # 26, 14, 26
kw = {
'vmin':hycom_true.iloc[:, 2:].values.min(),
'vmax':hycom_true.iloc[:, 2:].values.max(),
'levels':np.linspace(hycom_true.iloc[:, 2:].values.min(), hycom_true.iloc[:, 2:].values.max(), 10),
'alpha':0.7,
'cmap':'hot',
}
a = ax.contourf(X, Y, hycom_true[0].values.reshape(-1, 14), zdir='z', offset=0, **kw)
b = ax.contourf(X, Y, hycom_true[500].values.reshape(-1, 14), zdir='z', offset=-500, **kw)
# Set limits of the plot from coord limits
xmin, xmax = X.min(), X.max()
ymin, ymax = Y.min(), Y.max()
zmin, zmax = -z[-1], z[0]
ax.set(xlim=[xmin, xmax], ylim=[ymin, ymax], zlim=[zmin, zmax])
ax.set(
xlabel='Longitude/°E',
ylabel='Latitude/°N',
zlabel='Z/m',
zticks= [0, -100, -200, -500, -800],
zticklabels = [0, 100, 200, 500, 800])
cl = fig.colorbar(a, fraction=0.02, pad=0.1)
plt.show()
浙公网安备 33010602011771号