使用 Basemap、matplotlib 和轮廓命令绘制 netCDF 数据时出现 IndexError

Question

我正在尝试从 MERRA 再分析数据的 netcdf 文件在单个级别上绘制压力和风速变量。我正在使用底图模块进行绘图和其他必要的包来获取数据。不幸的是，我最终遇到了一个错误。

这是我的代码：

from netCDF4 import Dataset as NetCDFFile 
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.basemap import Basemap
import os
import netCDF4
from netCDF4 import Dataset
from osgeo import gdal 

os.chdir('F:\Atmospheric rivers\MERRA')
dirrec=str('F:\Atmospheric rivers\MERRA')

ncfile='MERRA2_400.inst6_3d_ana_Np.20180801.SUB.nc'
file=Dataset(ncfile,'r',format='NETCDF4')
lon=file.variables['lon'][:]
lat=file.variables['lat'][:]
time=file.variables['time'][:]

u=file.variables['U'][:]
v=file.variables['V'][:]
p=file.variables['PS'][:]
q=file.variables['QV'][:]

map = Basemap(projection='merc',llcrnrlon=70.,llcrnrlat=8.,urcrnrlon=90.,urcrnrlat=20.,resolution='i')

map.drawcoastlines()
map.drawstates()
map.drawcountries()
#map.drawlsmask(land_color='Linen', ocean_color='#CCFFFF', resolution ='i') # can use HTML names or codes for colors
#map.drawcounties() 
parallels = np.arange(0,50,5.) # make latitude lines ever 5 degrees from 30N-50N
meridians = np.arange(-95,-70,5.) # make longitude lines every 5 degrees from 95W to 70W
map.drawparallels(parallels,labels=[1,0,0,0],fontsize=10)
map.drawmeridians(meridians,labels=[0,0,0,1],fontsize=10)
x,y= np.meshgrid(lon-180,lat) # for this dataset, longitude is 0 through 360, so you need to subtract 180 to properly display on map
#x,y = map(lons,lats)
clevs = np.arange(960,1040,4)
cs = map.contour(x,y,p[0,:,:]/100,clevs,colors='blue',linewidths=1.)


plt.show()
plt.savefig('Pressure.png')

错误

%run "F:/Atmospheric rivers/MERRA/aosc.py"
C:\Users\Pavilion\AppData\Local\Enthought\Canopy\edm\envs\User\lib\site-packages\mpl_toolkits\basemap\__init__.py:3505: MatplotlibDeprecationWarning: The ishold function was deprecated in version 2.0.
  b = ax.ishold()
C:\Users\Pavilion\AppData\Local\Enthought\Canopy\edm\envs\User\lib\site-packages\mpl_toolkits\basemap\__init__.py:3570: MatplotlibDeprecationWarning: axes.hold is deprecated.
    See the API Changes document (http://matplotlib.org/api/api_changes.html)
    for more details.
  ax.hold(b)
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
F:\Atmospheric rivers\MERRA\aosc.py in <module>()
     36 #x,y = map(lons,lats)
     37 clevs = np.arange(960,1040,4)
---> 38 cs = map.contour(x,y,p[0,:,:]/100,clevs,colors='blue',linewidths=1.)
     39 
     40 
C:\Users\Pavilion\AppData\Local\Enthought\Canopy\edm\envs\User\lib\site-packages\mpl_toolkits\basemap\__init__.py in with_transform(self, x, y, data, *args, **kwargs)
    519             # convert lat/lon coords to map projection coords.
    520             x, y = self(x,y)
--> 521         return plotfunc(self,x,y,data,*args,**kwargs)
    522     return with_transform
    523 
C:\Users\Pavilion\AppData\Local\Enthought\Canopy\edm\envs\User\lib\site-packages\mpl_toolkits\basemap\__init__.py in contour(self, x, y, data, *args, **kwargs)
   3540                 # only do this check for global projections.
   3541                 if self.projection in _cylproj + _pseudocyl:
-> 3542                     xx = x[x.shape[0]/2,:]
   3543                     condition = (xx >= self.xmin) & (xx <= self.xmax)
   3544                     xl = xx.compress(condition).tolist()
IndexError: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices

所以，总而言之，错误结果是

IndexError：只有整数、切片 (:)、省略号 (...)、numpy.newaxis (None) 和整数或布尔数组是有效的索引

我该如何解决这个错误？

Answer 1

我猜你没有使用正确的 x 和 y 值来绘制绘图。底图轮廓要求所有 x、y 和 z 都是二维数组。在您的情况下， x 和 y 最有可能是经度和纬度值的向量。另外，需要根据地图坐标/投影，将经纬度转换为图形坐标。

你可以试试：

lons,lats = np.meshgrid(lon-180.,lat);
x,y = map(lons,lats)
cs = map.contour(x,y,p[0,:,:]/100,clevs,colors='blue',linewidths=1.)

第二行应将经度和纬度值转换为所选投影的正确 x 和 y 值。