0

我有一个模型级别的风速 netcdf。在同一个 netcdf 上，我有每个模型级别的高度。我把netcdf转换成一个立方体，所以每一层的高度就变成了一个辅助坐标。我想绘制一个横截面（经度 x 经度），并希望模型水平遵循地形。我尝试使用 Iris 模块文档示例（https://scitools.org.uk/iris/docs/latest/examples/General/cross_section.html），但它不起作用。

由于我已经有了每个级别相对于海平面的高度，我只需要切成垂直部分并绘制。我尝试按照文档示例进行操作，但它返回错误： ValueError: shape mismatch: objects cannot be broadcast to a single shape

贝娄是 Xarray.Dataset：

<xarray.Dataset>
Dimensions:    (latitude: 49, level: 21, longitude: 49)
Coordinates:
  * longitude  (longitude) float32 -52.0 -51.75 -51.5 ... -40.5 -40.25 -40.0
  * latitude   (latitude) float32 -15.0 -15.25 -15.5 ... -26.5 -26.75 -27.0
  * level      (level) int32 1 2 3 4 5 6 7 8 9 10 ... 13 14 15 16 17 18 19 20 21
    time       datetime64[ns] 2017-01-01T10:00:00
    altitude   (level, latitude, longitude) float32 271.3289 ... 1142.3843
Data variables:
    ws         (level, latitude, longitude) float32 0.77094275 ... 14.978188

我将 DataArray 转换为名为 ws_iris 的多维数据集：

Ws (unknown)    level   latitude    longitude
Shape            21       49              49
Dimension coordinates           
level             x       -              -
latitude        -         x              -
longitude       -         -              x
Auxiliary coordinates           
altitude       x          x              x
Scalar coordinates          
time    2017-01-01 10:00:00

这是我的代码：

import iris
import iris.plot as iplt
import iris.quickplot as qplt
import xarray as xr

ws = xr.open_dataset('ws.nc')

ws_iris = ws.ws.to_iris()

cross_section = next(ws_iris.slices(['longitude', 'level']))

qplt.contourf(cross_section, coords=['longitude', 'altitude'], cmap='viridis', levels=10)

如您所见，我遵循与文档示例相同的步骤（https://scitools.org.uk/iris/docs/latest/examples/General/cross_section.html）

当我绘制相对于模型级别的经度时，会绘制垂直部分：

qplt.contourf (cross_section, coords = ['longitude', 'level'], cmap = 'viridis', levels = 10)

但是当我绘制与高度相关的图时：

qplt.contourf (cross_section, coords = ['longitude', 'altitude'], cmap = 'viridis', levels = 10)

我收到以下错误消息：

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-35-c93b92f5c581> in <module>
----> 1 qplt.contourf(cross_section, coords=['longitude', 'altitude'], cmap='viridis', levels=10)

~\Anaconda3\lib\site-packages\iris\quickplot.py in contourf(cube, *args, **kwargs)
    202     coords = kwargs.get('coords')
    203     axes = kwargs.get('axes')
--> 204     result = iplt.contourf(cube, *args, **kwargs)
    205     _label_with_points(cube, result, coords=coords, axes=axes)
    206     return result

~\Anaconda3\lib\site-packages\iris\plot.py in contourf(cube, *args, **kwargs)
    937     coords = kwargs.get('coords')
    938     kwargs.setdefault('antialiased', True)
--> 939     result = _draw_2d_from_points('contourf', None, cube, *args, **kwargs)
    940 
    941     # Matplotlib produces visible seams between anti-aliased polygons.

~\Anaconda3\lib\site-packages\iris\plot.py in _draw_2d_from_points(draw_method_name, arg_func, cube, *args, **kwargs)
    516 
    517         u, v = plot_arrays
--> 518         u, v = _broadcast_2d(u, v)
    519 
    520         axes = kwargs.pop('axes', None)

~\Anaconda3\lib\site-packages\iris\plot.py in _broadcast_2d(u, v)
    234     u = np.atleast_2d(u)
    235     v = np.atleast_2d(v.T).T
--> 236     u, v = np.broadcast_arrays(u, v)
    237     return u, v
    238 

~\Anaconda3\lib\site-packages\numpy\lib\stride_tricks.py in broadcast_arrays(*args, **kwargs)
    257     args = [np.array(_m, copy=False, subok=subok) for _m in args]
    258 
--> 259     shape = _broadcast_shape(*args)
    260 
    261     if all(array.shape == shape for array in args):

~\Anaconda3\lib\site-packages\numpy\lib\stride_tricks.py in _broadcast_shape(*args)
    191     # use the old-iterator because np.nditer does not handle size 0 arrays
    192     # consistently
--> 193     b = np.broadcast(*args[:32])
    194     # unfortunately, it cannot handle 32 or more arguments directly
    195     for pos in range(32, len(args), 31):

ValueError: shape mismatch: objects cannot be broadcast to a single shape

有人可以帮我确定我的代码有什么问题吗？

或者也许有人知道如何制作我想要的图表。

非常感谢您提前。

4

1 回答 1

1

我对 iris 不熟悉，所以我想建议一种替代方法来绘制横截面xarray.Datatset。

Metpy 提供了一个非常有用的功能。对于这个问题。

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

否则，您可以使用以下代码对数据阵列进行切片：

    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))

和

def geodesic(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
于 2019-10-11T14:35:39.630 回答