Plotting 2d array on a sphere

[dennissergeev]

Description

Does pyvista allow for "stretching" a 2D array over a sphere?
(Related to #52, but simpler.)

I would like to create a 3D sphere with a 2D array as its surface. I tried creating a VTK sphere with specific resolution and then filling cells with flattened data, e.g.

import pyvista as pv

mesh = pv.Sphere(theta_resolution=45, phi_resolution=75, start_phi=1, end_phi=179)
mesh.cell_arrays["data"] = my_array.flatten(order="F")

However, I don't understand what the correct the number of cells / points depending on resolution should be, so my data dimensions don't match and in the end I end up with a mess. Is this because VTK sphere has tiangles as cells? What is the correct order of a data array then?

What I would like to get, can be shown using matplotlib and cartopy:

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import xarray as xr

lsm = xr.open_dataarray("lsm_4x5.nc")  # netCDF file is attached as a tar.gz archive

# Figure 1
lsm.plot()

# Figure 2
fig, ax = plt.subplots(
    subplot_kw=dict(projection=ccrs.Orthographic(central_longitude=180, central_latitude=45))
)
ax.imshow(lsm, transform=ccrs.PlateCarree())

Thanks!

Example Data

lsm_4x5.tar.gz

[banesullivan]

What you'll need to do is create a mesh and transform it to a sphere - the dataset you give isn't a complete sphere (latitudes range from -88 to 88 and longitudes range from 0 to 355), so you cannot inject his specific dataset onto a PolyData sphere object without significantly customizing it. Considering the dataset you give has latitudes and longitudes, I'd avoid creating a Sphere surface and just make a grid of the given spatial reference like so:

import xarray as xr
import numpy as np
import pyvista as pv

lsm = xr.open_dataarray("lsm_4x5.nc")  # netCDF file is attached as a tar.gz archive

xx, yy, zz = np.meshgrid(np.radians(lsm['longitude']), 
                         np.radians(lsm['latitude']), 
                         [0])

# Transform to spherical coordinates
radius = 6371.0e6
x = radius * np.cos(yy) * np.cos(xx)
y = radius * np.cos(yy) * np.sin(xx)
z = radius * np.sin(yy)

grid = pv.StructuredGrid(x, y, z)
grid['lsm'] = np.array(lsm).ravel(order='F')

grid.plot()

[banesullivan]

I suppose you might actually want this data to be on the cells and not the vertices of the mesh... in that case, maybe try this which will create a closed sphere:

(however this requires some knowledge about the spatial reference which could likely be implemented programatically besides explicitly setting the arange values)

import xarray as xr
import numpy as np
import pyvista as pv

lsm = xr.open_dataarray("lsm_4x5.nc")  # netCDF file is attached as a tar.gz archive

xx, yy, zz = np.meshgrid(np.radians(np.arange(0, 365, 5, )), 
                         np.radians(np.arange(-90, 94, 4)), 
                         [0])

# Transform to spherical coordinates
radius = 6371.0e6
x = radius * np.cos(yy) * np.cos(xx)
y = radius * np.cos(yy) * np.sin(xx)
z = radius * np.sin(yy)

grid = pv.StructuredGrid(x, y, z)
grid.cell_arrays['lsm'] = np.array(lsm).ravel(order='F')

grid.plot()

[dennissergeev]

This is brilliant, thank you!

[dennissergeev]

Continuing this topic, I was wondering if it is also possible to create a spherical "shell", i.e. a mesh that has multiple z levels instead of just [0]?

[banesullivan]

Hm... maybe it is possible to do it with some meshgrid tricks. I’ll see if I can make something work. I remember trying this back when I was looking at #52, but didn’t figure it out then

_{Sent with GitHawk}

[banesullivan]

Check this out and let me know what you think:

import numpy as np
import pyvista as pv

radial = np.linspace(0, 6371.0e6, 10)
xx, yy, zz = np.meshgrid(np.radians(np.arange(0, 365, 5, )), 
                         np.radians(np.arange(-90, 94, 4)), 
                         radial)

# Transform to spherical coordinates
x = zz * np.cos(yy) * np.cos(xx)
y = zz * np.cos(yy) * np.sin(xx)
z = zz * np.sin(yy)

grid = pv.StructuredGrid(x, y, z)
grid['radials'] = zz.ravel(order='F')

# Show wireframe to see inside of mesh
grid.wireframe().plot(notebook=0)

[dennissergeev]

Thanks a lot, @banesullivan! This is great.

This is exactly what I meant, so the next step that I tried is visualising something in this 3D wireframe, for example an isosurface. It seems to work for me, but only if I transpose the data array, i.e.:

grid.cell_arrays["iso_arr"] = scalar3d_data.transpose((1, 2, 0)).ravel(order="F")

where scalar3d has dimensions (height, latitude, longitude). The order of axes looks odd, but it seems to give what I expect. Does this make sense to you?

[banesullivan]

That order has to do with the F ordering (z, y, x) so when you ravel it it comes out right

[banesullivan]

or maybe not, I'm not sure. This is a bit messy

[dennissergeev]

OK, thanks!