Я использую Cartopy для построения данных модели над Антарктидой (файл netCDF — ссылка для загрузки приведена ниже). В этом случае я строю контуры высотой 500 мб. Проблема: данные не являются непрерывными на 180 градусов на восток, и, поскольку я рисую в полярной стереографической проекции, мне нужно, чтобы долгота была непрерывной. На рисунке вы можете видеть, что контуры отрисовываются вокруг круга широты на 180 градусов долготы с другой стороны вместо того, чтобы соединяться поперек на 180 градусах.
Я нашел вспомогательную функцию Cartopy add_cyclo_point cartopy.util.add_cyclo_point(data, coord=None, axis=-1), чтобы попытаться сделать долготу непрерывной (и мне удалось использовать это решение в другом наборе данных); однако это не работает с этим набором данных по двум причинам, из которых я могу сделать вывод:
- Данные широты и долготы находятся в двумерной сетке, а координата должна быть одномерным массивом.
- Значения координат (в моем случае долготы) должны быть расположены через равные промежутки.
Для 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()