Как построить вертикальный разрез атмосферы вместе с топографией с помощью модуля Iris? - PullRequest
Новая
       64

Как построить вертикальный разрез атмосферы вместе с топографией с помощью модуля Iris?

0 голосов
/ 10 октября 2019

У меня есть netcdf со скоростью ветра на модельных уровнях. В том же netcdf у меня есть высота каждого уровня модели. Я преобразовал netcdf в куб, поэтому высота каждого уровня стала вспомогательной координатой. Я хотел бы построить поперечное сечение (долгота х долгота) и хотел бы, чтобы уровни модели соответствовали топографии. Я попытался использовать пример документации модуля Iris (https://scitools.org.uk/iris/docs/latest/examples/General/cross_section.html),, но он не работает.

Поскольку у меня уже есть высота каждого уровня относительно уровня моря, мне просто нужно разрезать на вертикальные участки и построить график. Я попытался следовать примеру документации, но он возвращает ошибку: ValueError: несоответствие формы: объекты не могут быть переданы одной фигуре

Ниже - 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)

Iполучите следующее сообщение об ошибке: 

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

Может кто-нибудь помочь мне определить, что не так с моим кодом.

Или, может быть, кто-то знает, как сделать график, который я хочу в противном случае.

спасибоочень заранее.

1 Ответ

0 голосов
/ 11 октября 2019

Я не знаком с радужной оболочкой, поэтому хочу предложить альтернативный способ построения поперечных сечений с помощью xarray.Datatset.

Metpy предоставляет очень полезную функцию. для этого вопроса.

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

В противном случае вы можете использовать следующий код для нарезки вашего 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))

С 

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