python - 使用 pyvista 和 pyproj 重新投影 vtu 文件的网格点
问题描述
有没有机会我可以用 ndarray 替换pyvista_ndarray的坐标?
我的目标是重新投影 vtu 文件(非结构化网格)的点坐标。网格对象的当前坐标位于坐标参考系统 EPSG 27672 中,我希望它们位于 EPSG 4326 (WGS84) 中。为此,我使用pyvista模块打开 vtu 文件:
mesh = pv.read("arg_Tm2__t0002.pvtu")
type(mesh.points)
>pyvista.core.pyvista_ndarray.pyvista_ndarray
结果,mesh.points给出了 3 个空间坐标。然后,我使用pyproj模块将 3 个坐标重新投影到 EPSG 4326 中。通过组合生成的 3 个 x,y,z numpy.ndarray,我现在得到一个形状和大小类似于mesh.points的 NumPy 数组。
# set the pyproj transformer
crs1 = CRS.from_epsg(27562)
crs2 = CRS.from_epsg(4326)
reproj = Transformer.from_crs(crs1, crs2)
# Reprojection
reproj_dataY,reproj_dataX,altitude = reproj.transform(mesh.points[:,0],mesh.points[:,1],mesh.points[:,2])
reprojData = np.column_stack((reproj_dataX,reproj_dataY,altitude))
#compare objects
print('original Mesh points -> ', mesh.points)
print('original Mesh type: ', type(mesh.points))
print('Reprojected points-> ', reprojData)
print('Reprojected type: ', type(reprojData))
Original Mesh points -> [[958427.33 119680.95 2396.288549 ]
[957754.39 120023.85 2430.1833881 ]
[957256.56 120241.02 2112.22953263]
...
[963366.748527 115096.364632 3054.75408138]
[963401.840285 113351.753238 3024.50286566]
[963497.913738 113339.696062 3048.83674197]]
Original Mesh type: <class 'pyvista.core.pyvista_ndarray.pyvista_ndarray'>
Reprojected points-> [[ 6.96487903 45.9823843 2396.288549 ]
[ 6.95646936 45.98581994 2430.1833881 ]
[ 6.95021969 45.98803333 2112.22953263]
...
[ 7.02498443 45.93857775 3054.75408138]
[ 7.02409542 45.92289079 3024.50286566]
[ 7.02532248 45.92273099 3048.83674197]]
Reprojected type: <class 'numpy.ndarray'>enter code here
现在,是时候替换 vtu 对象的坐标了:
mesh.points = reprojData
最后,我检查了修改后的网格:X 边界和Y 边界已被修改并且范围是正确的。但是,该图显示的是一行点,而不是一个漂亮的 3d 对象。:(。
你知道有什么问题吗?您是否看到另一种管理重投影的方法?
解决方案
变换后 XY 和 Z 的取值范围有显着差异:
>>>np.array(mesh.bounds).reshape((3,-1)).ptp(axis=1)
array([1.22515302e-01, 7.78657599e-02, 2.47978788e+03])
XY 确实以度为单位,Z 仍然以米为单位。这些数据的视觉表示是不切实际的,因此应该调整数据。
推荐阅读
- r - 计算子组的非 NA 值的百分比
- ruby-on-rails - 升级到 Rails 6 后参数数量错误(给定 4,预期为 0..1)
- mongodb - 使用 MongoDB 处理时间间隔的最佳方法是什么
- postgresql - RDS Postgresql CPU 使用率 100%,未知连接“其他”
- javascript - 如何使用 Javascript 从 hls(m3u8) 获取 mp3 文件?
- php - 正则表达式 该帖子首先出现在
- javascript - 如何在花括号内的表达式周围放置强标签?
- database - Sqlite中B-Tree的度数是多少?
- python - 将同一行中的多个值连接到列表中
- c++ - OpenMP 和 GMP 瓶颈问题