How To Plot A Vertical Section Of The Atmosphere Along With The Topography Using The Iris Module?

I have a netcdf with wind speed at model levels. On the same netcdf I have the altitude of each model level. I converted netcdf to a cube, so the altitude of each level became an a

Solution 1:

I am not familiar with iris so I want to suggest an alternative way to plot cross sections with xarray.Datatset.

Metpy provides a very helpful function. for this issue.

https://unidata.github.io/MetPy/latest/api/generated/metpy.interpolate.cross_section.html

Otherwise you can use following code to slice your datarray:

    points_cross = geodesic(crs_data, start, end, steps)

    data_sliced = data.interp({
        x.name: xr.DataArray(points[:, 0], dims='index', attrs=x.attrs),
        y.name: xr.DataArray(points[:, 1], dims='index', attrs=y.attrs)
    }, method='linear')
    data_sliced.coords['index'] = range(len(points))

With

defgeodesic(crs, start, end, steps):
    r"""Construct a geodesic path between two points.

    This function acts as a wrapper for the geodesic construction available in `pyproj`.

    Parameters
    ----------
    crs: `cartopy.crs`
        Cartopy Coordinate Reference System to use for the output
    start: (2, ) array_like
        A latitude-longitude pair designating the start point of the geodesic (units are
        degrees north and degrees east).
    end: (2, ) array_like
        A latitude-longitude pair designating the end point of the geodesic (units are degrees
        north and degrees east).
    steps: int, optional
        The number of points along the geodesic between the start and the end point
        (including the end points).

    Returns
    -------
    `numpy.ndarray`
        The list of x, y points in the given CRS of length `steps` along the geodesic.

    See Also
    --------
    cross_section

    """import cartopy.crs as ccrs
    from pyproj import Geod

    g = Geod(crs.proj4_init)
    geodesic = np.concatenate([
        np.array(start[::-1])[None],
        np.array(g.npts(start[1], start[0], end[1], end[0], steps - 2)),
        np.array(end[::-1])[None]
    ]).transpose()
    points = crs.transform_points(ccrs.Geodetic(), *geodesic)[:, :2]

    return points