Plots With PyVista

From Planets
Revision as of 09:19, 4 September 2022 by Martin.turbet (talk | contribs)

Jump to: navigation, search

We provide here a tutorial on how to make pretty visuals from GCM simulations. Do not hesitate to send us your best visuals so that we can enrich our GCM plot gallery.


import pyvista as pv
import numpy as np
from netCDF4 import Dataset
# import all the required libraries (netcdf4, pyvista)
#FIRST WE GET THE DATA FROM THE GCM SIMULATION

# 4D GRID VARIABLES
nc1 = Dataset('diagfi.nc')

Time_GCM=nc1.variables['Time'][:]
lat_GCM=nc1.variables['latitude'][:]
lon_GCM=nc1.variables['longitude'][:]
alt_GCM=nc1.variables['altitude'][:]
aire_GCM=nc1.variables['aire'][:][:]

resolution_lon_GCM = abs(lon_GCM[1]-lon_GCM[0])
resolution_lat_GCM = abs(lat_GCM[1]-lat_GCM[0])

# 2D VARIABLES
tsurf1=nc1.variables['tsurf'][:][:][:]
OLR1=nc1.variables['OLR'][:][:][:]
h2o_ice_col1=nc1.variables['h2o_ice_col'][:][:][:]

# 3D VARIABLES
u1=nc1.variables['u'][:][:][:][:]
v1=nc1.variables['v'][:][:][:][:]
w1=nc1.variables['w'][:][:][:][:]
ice1=nc1.variables['h2o_ice'][:][:][:][:]
temp1=nc1.variables['temp'][:][:][:][:]
# DEFINE REFERENCE RADIUS OF THE PLANET
RADIUS = 1.0

# Longitudes and latitudes used to plot the GCM simulation
x = lon_GCM[0:len(lon_GCM)-1]+180.
y = lat_GCM[0:len(lat_GCM)]

# Define polar coordinates
y_polar = 90.0 - y  # grid_from_sph_coords() expects polar angle


# FUNCTION TO CALCULATE COORDINATE CELL BOUNDARIES
def _cell_bounds(points, bound_position=0.5):
    """
    Calculate coordinate cell boundaries.

    Parameters
    ----------
    points: numpy.array
        One-dimensional array of uniformly spaced values of shape (M,)
    bound_position: bool, optional
        The desired position of the bounds relative to the position
        of the points.

    Returns
    -------
    bounds: numpy.array
        Array of shape (M+1,)

    Examples
    --------
    >>> a = np.arange(-1, 2.5, 0.5)
    >>> a
    array([-1. , -0.5,  0. ,  0.5,  1. ,  1.5,  2. ])
    >>> cell_bounds(a)
    array([-1.25, -0.75, -0.25,  0.25,  0.75,  1.25,  1.75,  2.25])
    """
    assert points.ndim == 1, "Only 1D points are allowed"
    diffs = np.diff(points)
    delta = diffs[0] * bound_position
    bounds = np.concatenate([[points[0] - delta], points + delta])
    return bounds


# Create arrays of grid cell boundaries, which have shape of (x.shape[0] + 1)
xx_bounds = _cell_bounds(x)
yy_bounds = _cell_bounds(y_polar)