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 |
| |
| |
| 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) |
| |
| |
| lons = hycom['lon'].values |
| lats = hycom['lat'].values |
| depth = hycom['depth'].values |
| |
| |
| 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) |
注意上述程序因为插值后的网格点太多，我没跑出来结果，应该是时间不够，因为我把插值后的网格改少之后是可以跑通的。比如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) |
| 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) |
这个速度很快，不过效果略差，感觉不如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) |
| |
| |
| lons = hycom_['lon'].values |
| lats = hycom_['lat'].values |
| temp = hycom_[0].values |
| |
| |
| OK = OrdinaryKriging(lons, lats, temp, variogram_model='linear') |
| z1, ss1 = OK.execute('grid', grid_lon, grid_lat) |
| print(z1.shape) |
插值速度比cubic慢一些，不过也挺快的了，效果也很好，不过和cubic（三次样条插值-二维）肉眼看差距不大
1.2.2.scipy的griddata的cubic（三次样条插值）
| from scipy.interpolate import griddata |
| lat_lon_points = np.array(hycom_[['lat', 'lon']], dtype=np.float32) |
| grid_x, grid_y = np.mgrid[20.53:21.52:100j, 119.01:120:100j] |
| values = data[0].values |
| values_inter = griddata(lat_lon_points, values, (grid_x, grid_y), method='cubic') |
| print(values_inter.shape) |
效果看着和克里格没啥区别，速度也是
1.2.3.scipy的griddata的linear二维插值
| from scipy.interpolate import griddata |
| lat_lon_points = np.array(hycom_[['lat', 'lon']], dtype=np.float32) |
| grid_x, grid_y = np.mgrid[20.53:21.52:100j, 119.01:120:100j] |
| values = data[0].values |
| values_inter = griddata(lat_lon_points, values, (grid_x, grid_y), method='linear') |
| print(values_inter.shape) |
效果看起来也可以，不过感觉比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) |
效果肉眼上介于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) |
| |
| 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) |
| 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) |
| |
| 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) |
| 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_xlabel('Longitude/°E') |
| ax.set_ylabel('Latitude/°N') |
| |
| ax.grid(which='both', axis='both', alpha=0.4) |
| |
| |
| plt.legend(loc='upper right') |
| |
| |
| 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] |
| 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() |
| 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)] |
| |
| |
| ) |
| |
| 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) |
| |
| 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) |
| |
| |
| 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() |
2.2.3.三维等温线图contourf（）展示海水在垂直方向上温度变化趋势，就是那种分层展示的，不过感觉用处不大
2.2.4.三维立体像素图，voxels（），用来显示海水在垂直方向上温度变化趋势，以及每层的变化趋势
【推荐】国内首个AI IDE，深度理解中文开发场景，立即下载体验Trae
【推荐】编程新体验，更懂你的AI，立即体验豆包MarsCode编程助手
【推荐】抖音旗下AI助手豆包，你的智能百科全书，全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell：AI 加持，快人一步
· 25岁的心里话
· 闲置电脑爆改个人服务器（超详细） #公网映射 #Vmware虚拟网络编辑器
· 零经验选手，Compose 一天开发一款小游戏！
· 通过 API 将Deepseek响应流式内容输出到前端
· 因为Apifox不支持离线，我果断选择了Apipost！