0

Is there a chance I can replace coordinates of a pyvista_ndarray by a ndarray?

My objective is to reproject points coordinates of a vtu file (unstructured grid). The current coordinates of the mesh object are in the coordinate reference system EPSG 27672 and I want them in EPSG 4326 (WGS84). For this, I open the vtu file using the pyvista module :

mesh = pv.read("arg_Tm2__t0002.pvtu")
type(mesh.points)
>pyvista.core.pyvista_ndarray.pyvista_ndarray

As a result, mesh.points gives the 3 spatial coordinates. Then, I use pyproj module to reproject the 3 coordinates into EPSG 4326. By combining the 3 resulting x,y,z numpy.ndarray, I now get a NumPy array with a shape and size similar to mesh.points.

# 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

Now, It's time to replace the coordinates of the vtu object:

mesh.points = reprojData

Finally, I check the modified mesh: The X bounds and Y bounds have been modified and the ranges are correct. However, the plot shows a line of points instead a nice 3d object. :(.

Do you have any idea what is wrong? Do you see another way to manage reprojection?

jean Baba
  • 11
  • 4
  • `mesh.points = reprojData` should already work. `points` is a property that generates a new `pyvista_ndarray` each time, and assigning an array of shape `(mesh.n_points, 3)` to `mesh.points` should indeed overwrite the point coordinates with the new ones. Are you sure the problem is not coming from elsewhere? And what is "the plot" that shows a line? – Andras Deak -- Слава Україні Apr 06 '21 at 20:22
  • 1
    thanks @AndrasDeak. As you said, the problem is coming from elsewhere : the difference in units after the reprojection: X and Y are in degrees and Z still in meter. This leads to a very large difference in the range of values between XY and Z. As a result, the plot is unfeasible. – jean Baba Apr 07 '21 at 07:47

1 Answers1

1

The range of values of XY and Z after the transformation are significantly different:

>>>np.array(mesh.bounds).reshape((3,-1)).ptp(axis=1)
array([1.22515302e-01, 7.78657599e-02, 2.47978788e+03])

XY are indeed in degrees and Z still in meter. The visual representation of this data is unrealistic and data should be adapted to do so.

jean Baba
  • 11
  • 4