Использование add_cyclo_point в Cartopy: ошибка: координаты должны быть расположены через равные промежутки

Я использую Cartopy для построения данных модели над Антарктидой (файл netCDF — ссылка для загрузки приведена ниже). В этом случае я строю контуры высотой 500 мб. Проблема: данные не являются непрерывными на 180 градусов на восток, и, поскольку я рисую в полярной стереографической проекции, мне нужно, чтобы долгота была непрерывной. На рисунке вы можете видеть, что контуры отрисовываются вокруг круга широты на 180 градусов долготы с другой стороны вместо того, чтобы соединяться поперек на 180 градусах.

Я нашел вспомогательную функцию Cartopy add_cyclo_point cartopy.util.add_cyclo_point(data, coord=None, axis=-1), чтобы попытаться сделать долготу непрерывной (и мне удалось использовать это решение в другом наборе данных); однако это не работает с этим набором данных по двум причинам, из которых я могу сделать вывод:

  1. Данные широты и долготы находятся в двумерной сетке, а координата должна быть одномерным массивом.
  2. Значения координат (в моем случае долготы) должны быть расположены через равные промежутки.

Для 1: я думаю, что могу решить эту проблему, просто сгладив 2D-массив? Для 2: я нашел эту ссылку, чтобы в целом объяснить проблему, с которой я столкнулся. вероятно, имеет, но, поскольку мои данные о долготе в 2D, я не уверен, что точно следую объяснению.

Вот фрагмент кода, в котором я пытаюсь использовать функцию add_cyclo_point:

import netCDF4 as nc 
import numpy as np
from cartopy.util import add_cyclic_point

file = '2016061000_WRF_d2_HGT_500_f006.nc'
nc_fid  = nc.Dataset(file,'r')

height = np.array(nc_fid.variables['HGT'])
lons= np.array(nc_fid.variables['g5_lon_1'])

lons_flat = lons.flatten()
height_flat = height.flatten()

heightC, lonC = add_cyclic_point(height_flat, coord=lons_flat)

Вот ошибка, которую он выдает:


---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-37-214135850083> in <module>
----> 1 heightC, lonC = add_cyclic_point(height_flat, coord=lons_flat)

~/miniconda3/lib/python3.7/site-packages/cartopy/util.py in add_cyclic_point(data, coord, axis)
     99         delta_coord = np.diff(coord)
    100         if not np.allclose(delta_coord, delta_coord[0]):
--> 101             raise ValueError('The coordinate must be equally spaced.')
    102         new_coord = ma.concatenate((coord, coord[-1:] + delta_coord[0]))
    103     slicer = [slice(None)] * data.ndim

ValueError: The coordinate must be equally spaced.

Вот метаданные из файла:

OrderedDict([('g5_lon_1', <class 'netCDF4._netCDF4.Variable'>
              float32 g5_lon_1(y, x)
                  standard_name: longitude
                  long_name: longitude
                  units: degrees_east
                  _CoordinateAxisType: Lon
              unlimited dimensions: 
              current shape = (627, 666)
              filling on, default _FillValue of 9.969209968386869e+36 used),
             ('g5_lat_0', <class 'netCDF4._netCDF4.Variable'>
              float32 g5_lat_0(y, x)
                  standard_name: latitude
                  long_name: latitude
                  units: degrees_north
                  _CoordinateAxisType: Lat
              unlimited dimensions: 
              current shape = (627, 666)
              filling on, default _FillValue of 9.969209968386869e+36 used),
             ('HGT', <class 'netCDF4._netCDF4.Variable'>
              float32 HGT(y, x)
                  long_name: Geopotential height
                  units: gpm
                  coordinates: g5_lon_1 g5_lat_0
                  _FillValue: 1e+20
                  missing_value: 1e+20
                  initial_time: 06/10/2016 (00:00)
                  forecast_time_units: hours
                  forecast_time: 6
                  level: 500
                  parameter_number: 7
                  parameter_table_version: 2
                  gds_grid_type: 5
                  level_indicator: 100
                  center: National Center for Atmospheric Research (NCAR), Boulder, CO
              unlimited dimensions: 
              current shape = (627, 666)
              filling on)])

Я предполагаю, что это связано с сеткой, в которой находятся данные (учитывая, что данные широты и долготы находятся в двумерной сетке). В моем другом наборе данных были одномерные данные широты и долготы.

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

file = '2016061000_WRF_d2_HGT_500_f006.nc'
nc_fid  = nc.Dataset(file,'r')

height = nc_fid.variables['HGT']
lons= nc_fid.variables['g5_lon_1']
lats = nc_fid.variables['g5_lat_0']


fig = plt.figure(figsize=(5,5))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.SouthPolarStereo())
land_res = cfeature.NaturalEarthFeature('physical','land','10m',edgecolor='face',facecolor=cfeature.COLORS['land'])

ax.contour(lons[:,:],lats[:,:],height[:,:],linewidths=0.5,colors='black',transform=ccrs.PlateCarree())

ax.add_feature(cfeature.COASTLINE.with_scale('10m'),linewidth=2.0)
ax.add_feature(cfeature.BORDERS)
ax.add_feature(cfeature.STATES)
ax.add_feature(land_res)
ax.gridlines(alpha=.5)
plt.show()

netCDFfile

Карта Cartopy с контурами высотой 500 МБ, построенными на основе данных. Я предоставил файл netCDF по предоставленной ссылке.


person wxstudent    schedule 07.08.2020    source источник