conversion from pv.mesh to gv.mesh?

[k-a-mendoza]

📚 Documentation

How does one add a mesh made in normal pyvista to a geovista plotter?

for example, I have a draped 2D line id like to place over a pv sphere. Its already in x,y lat lon coords, but it also has a z coord.

[welcome]

📢 Nice one! Your first issue! Thanks for telling us how to improve geovista 📢

[k-a-mendoza]

For those who want a quick and dirty way to do it and assuming your mesh is formatted as points = [lon, lat, elevation (km)], this is what I'm doing:

def convert_latlondepth_to_unit_circle(mesh: pv.DataSet):
    """converts a lat, lon, depth mesh to a unit circle"""
    # Get the min and max values of the mesh
    earth_radius = 6371.0
    lon = np.deg2rad(mesh.points[:,0])
    lat = np.deg2rad(mesh.points[:,1])
    elevation = 1 + (mesh.points[:,2] / earth_radius)
    # Convert spherical coordinates to cartesian coordinates
    x = elevation * np.cos(lat) * np.cos(lon)
    y = elevation * np.cos(lat) * np.sin(lon) 
    z = elevation * np.sin(lat)

    mesh.points = np.column_stack([x, y, z])
    return mesh

[bjlittle]

Hey @k-a-mendoza,

Great to hear from you 😄

I'm guessing that you have a workflow where you create the pyvista.DataSet instance already and you then want to cast it to another pyvista data structure e.g., pyvista.PolyData, but the points are serialized as lons, lats and elevation.

Perhaps you could try using the convenience of the geovista.common.to_cartesian function ...

e.g., something along the lines of:

from geovista.common import to_cartesian

def convert_latlondepth_to_unit_circle(mesh: pv.DataSet) -> pv.DataSet:
    radius = 6371.0
    zscale = 1 / radius
    lons = mesh.points[:, 0]
    lats = mesh.points[:, 1]
    elevation = mesh.points[:, 2]    

    mesh.points = to_cartesian(lons, lats, zlevel=elevation, zscale=zscale, stacked=True)

    return mesh

Or perhaps the geovista.bridge.Transform.from_points method might be closer to what you need?

from geovsta.bridge import Transform

lons = mesh.points[:, 0]
lats = mesh.points[:, 1]
elevation = mesh.points[0, 2]
zscale = 1 / 6371.0

mesh = Transform.from_points(lons, lats, zlevel=elevation, zscale=zscale)

The above will return a pyvista.PolyData point-cloud.

Does your pyvista.DataSet have any form of connectivity?

Let me know if this helps 👍

[k-a-mendoza]

@bjlittle that helps a lot! Although, it might be worth putting some functionality in the future to change the points directly without having to initially extract them and feed them in the right args and kwargs.

I have a few datasets actually; some point clouds as well as a logically-rectilinear mesh (although in the future I may add some unstructured meshes). The meshes usually come as a custom AEQD projection, so my process for adding them involves 1st using proj to convert them to latlons, then a cartesian conversion to get them into xyz format.

[bjlittle]

@k-a-mendoza Great, glad to help 🎉

However, I agree, my suggestions above are far from ideal. It's the first time I've seen lon/lat/elevation spatial data serialized directly within the points of a pyvista.DataSet.

My main focus up until now has been to mature the capability within the geovista.bridge as the main entry-point for users to create meshes from their various geospatial data.

Would it be possible to share with me some of your source datasets? How are they encoded? In NetCDF, GRIB or some other file format? Do you have example code of how you do this in a GH gist or repo? So that I can get a fuller understanding of what you're doing.

I'd be keen to make your workflow in geovista slicker and more natural, if possible 👍