-
-
Save dennissergeev/4cb8a4773f6d08c88ce4241dfe568236 to your computer and use it in GitHub Desktop.
LFRic wind vectors in GeoVista
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
import iris | |
from iris.experimental.ugrid import PARSE_UGRID_ON_LOAD | |
import numpy as np | |
import geovista as gv | |
import pyvista as pv | |
def cube2d_to_mesh(cube2d: iris.cube.Cube, **kwargs) -> pv.PolyData: | |
"""Construct a spherical mesh from a 2D iris cube using geovista.""" | |
mesh = gv.Transform.from_unstructured( | |
*[i.points for i in cube2d.mesh.node_coords], | |
cube2d.mesh.face_node_connectivity.indices_by_location(), | |
data=cube2d.data, | |
start_index=cube2d.mesh.face_node_connectivity.start_index, | |
**kwargs | |
) | |
return mesh | |
def transform_vectors_sph_to_cart(theta, phi, r, u, v, w): | |
"""Transform vectors from spherical (r, phi, theta) to cartesian coordinates (z, y, x).""" | |
# xx, yy, _ = np.meshgrid(np.radians(theta), np.radians(phi), r, indexing="ij") | |
# th, ph = xx.squeeze(), yy.squeeze() | |
sin_ph = np.sin(np.radians(phi)) | |
cos_ph = np.cos(np.radians(phi)) | |
sin_th = np.sin(np.radians(theta)) | |
cos_th = np.cos(np.radians(theta)) | |
ph = np.radians(phi) | |
# Transform wind components from spherical to cartesian coordinates | |
# https://en.wikipedia.org/wiki/Vector_fields_in_cylindrical_and_spherical_coordinates | |
u_t = sin_ph * cos_th * w + cos_ph * cos_th * v - sin_th * u | |
v_t = sin_ph * sin_th * w + cos_ph * sin_th * v + cos_th * u | |
w_t = cos_ph * w - sin_ph * v | |
return u_t, v_t, w_t | |
def main(): | |
with PARSE_UGRID_ON_LOAD.context(): | |
u, v, w = iris.load("thai_hab1_uvw_500d.nc") | |
u_t, v_t, w_t = transform_vectors_sph_to_cart( | |
u.coord("longitude").points, | |
u.coord("latitude").points, | |
gv.common.RADIUS, | |
u.data, | |
v.data, | |
w.data, | |
) | |
mesh = cube2d_to_mesh(u) | |
mesh["vectors"] = np.vstack((u_t, v_t, w_t)).T | |
mesh.set_active_vectors("vectors") | |
# mesh.arrows # <-- ValueError: data length of (13824) != required length (13826) | |
if __name__ == "__main__": | |
main() |
This file has been truncated, but you can view the full file.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
�HDF | |
�������� 0 n!N�OHDR�" |