使用loop从NetCDF生成多个底图

Question

我需要生成多个位势高度异常图（数据取自NCEP NCAR再分析），为此我使用python netCDF4和basemap包。我编写的代码生成单个图像，但希望对其进行修改以生成多个绘图。生成单个图像的代码如下所示：

from netCDF4 import Dataset, num2date
import os
import numpy as np
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
C=Dataset('G:\\GPANC\\hgt.2010.nc','r')
C3=C.variables['hgt'][177,2,:,:]
C1=C.variables['time']
dates = num2date(C1[:], C1.units)
dates[177:180]
X1 = [1970,1971,1973,1974,1975,1977,1978,1979,1981,1982,1983,1985,1986,1987,1989,1990,1991,1993,1994,1995,1997,1998,1999,2001,2002,2003,2005,2006,2007,2009,2010,2011,2013,2014,2015,2017,2018]
X2 = [1972,1976,1980,1984,1988,1992,1996,2000,2004,2008,2012,2016]
Z1=np.zeros([37,73,144])
row11=0
Z2=np.zeros([12,73,144])
row12=0
for year in X1:
  N1=Dataset('G:\\GPANC\\hgt.'+str(year)+'.nc','r')
  Z1[row11,:,:]=N1.variables['hgt'][177,2,:,:]
  row11=row11+1
  print(year)
for year in X2:
  N2=Dataset('G:\\GPANC\\hgt.'+str(year)+'.nc','r')
  Z2[row12,:,:]=N2.variables['hgt'][178,2,:,:]
  row12=row12+1
  print(year)
y=np.concatenate((Z1,Z2))
b1=N1.variables['lon'][:]
c1=N1.variables['lat'][:]
z=np.mean(y, axis=0)
S=C3-z
plt.subplot(111)
map = Basemap(projection='cyl',llcrnrlat= 4,urcrnrlat= 50,\
                 resolution='h', llcrnrlon=45,urcrnrlon=125)
map.drawparallels(np.arange(4,50,6),labels=[1,0,0,0],linewidth=0,fontsize=8)
map.drawmeridians(np.arange(45,125,8),labels=[0,0,0,1],linewidth=0,fontsize=8)
map.drawcountries(linewidth=1.5)
map.drawcoastlines(linewidth=1.5)
df=[{'lon': 80.3319, 'lat': 26.4499,'site':'Kanpur'}]
for point in df:
    u,v=map(point['lon'],point['lat'])
    plt.plot(u,v,marker='o',color='red',ms=5)
    plt.annotate(point['site'], xy=(u,v),  xycoords='data',xytext=(-8,3), textcoords = 'offset points')
x6,y6=np.meshgrid(b1,c1)
q6,p6=map(x6,y6)
CI=[-80,-72,-64,-56,-48,-40,-32,-24,-16,-8,0,8,16,24,32,40,48,56,64]
map.contourf(q6,p6,S,CI,cmap='gist_ncar',extend='both',zorder=1)
clb=plt.colorbar()
clb.set_label('m', labelpad=-40, y=1.05, rotation=0)
clb.set_ticks([-80,-72,-64,-56,-48,-40,-32,-24,-16,-8,0,8,16,24,32,40,48,56,64])
Z6=plt.contour(q6,p6,S,CI,colors='black',linewidths=1,zorder=2)
plt.clabel(Z6, inline=1, fontsize=8,fmt="%i")
plt.title('Geopotential height Anomalies (850 hPa) on 27-06-2010 (1970-2018 Climatology)',fontsize=10,pad=10)

#in order to develop multiple plots I made following changes
data=[]
for i in np.arange(177,180,1):
                C3=C.variables['hgt'][i,2,:,:]
                data.append(C3)
X1 = [1970,1971,1973,1974,1975,1977,1978,1979,1981,1982,1983,1985,1986,1987,1989,1990,1991,1993,1994,1995,1997,1998,1999,2001,2002,2003,2005,2006,2007,2009,2010,2011,2013,2014,2015,2017,2018]
X2 = [1972,1976,1980,1984,1988,1992,1996,2000,2004,2008,2012,2016]
Z1=np.zeros([37,73,144])
row11=0
Z2=np.zeros([12,73,144])
row12=0
for year in X1: 
               for i in np.arange(177,180,1):
                      N1=Dataset('G:\\GPANC\\hgt.'+str(year)+'.nc','r')
                      Z1[row11,:,:]=N1.variables['hgt'][i,2,:,:]
                      row11=row11+1
                      print(year)    # I get an error here
1970
1970
1970
1971
1971
1971
1973
1973
1973
1974
1974
1974
1975
1975
1975
1977
1977
1977
1978
1978
1978
1979
1979
1979
1981
1981
1981
1982
1982
1982
1983
1983
1983
1985
1985
1985
1986
Traceback (most recent call last):
  File "<stdin>", line 4, in <module>
IndexError: index 37 is out of bounds for axis 0 with size 37 

我这里有两个问题 1.如何存储多个NetCDF文件中日期177到180的数据，以及 2.如何使用存储的数据生成多个底图图，即177到180的所有日期[日期177表示2010年6月27日，以此类推]

Answer 1

谢谢你提供必要的链接，我能够生成情节。生成多个地质位高度异常图的代码如下所示：

#loading required packages
from netCDF4 import Dataset, num2date
import os
import numpy as np
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

C=Dataset('G:\\GPANC\\hgt.1977.nc','r')
data=[]
for i in np.arange(288,291,1):
                 C3=C.variables['hgt'][i,2,:,:]
                 data.append(C3)

X1 = [1970,1971,1973,1974,1975,1977,1978,1979,1981,1982,1983,1985,1986,1987,1989,1990,1991,1993,1994,1995,1997,1998,1999,2001,2002,2003,2005,2006,2007,2009,2010,2011,2013,2014,2015,2017,2018]
X2 = [1972,1976,1980,1984,1988,1992,1996,2000,2004,2008,2012,2016]

Z1=[]
for year in X1: 
            N1=Dataset('G:\\GPANC\\hgt.'+str(year)+'.nc','r')
            for i in np.arange(288,291,1):
                                    K1=N1.variables['hgt'][i,2,:,:]
                                    Z1.append(K1)
                                    print(year) 
Z2=[]
for year in X2: 
            N2=Dataset('G:\\GPANC\\hgt.'+str(year)+'.nc','r')
            for i in np.arange(289,292,1):
                                    K2=N2.variables['hgt'][i,2,:,:]
                                    Z2.append(K2)
                                    print(year) 
#non-leap year
P1 = [Z1[i] for i in np.arange(0,111,3)]
P2 = [Z1[i] for i in np.arange(1,111,3)]
P3 = [Z1[i] for i in np.arange(2,111,3)]

#leap year
Q1 = [Z2[i] for i in np.arange(0,36,3)]
Q2 = [Z2[i] for i in np.arange(1,36,3)]
Q3 = [Z2[i] for i in np.arange(2,36,3)]

#Concatenating the leap and non-leap year data
Y1=np.concatenate((P1,Q1))
Y2=np.concatenate((P2,Q2))
Y3=np.concatenate((P3,Q3))

#Anomaly
ANA1=data[0]-np.mean(Y1, axis=0)
ANA2=data[1]-np.mean(Y2, axis=0)
ANA3=data[2]-np.mean(Y3, axis=0)

#plotting with Basemap
#Coordinates of Thiruvananthapuram
lat = 8.5241
lon =76.9366
def plot_background(ax):
                  map = Basemap(projection='cyl',ax=ax,llcrnrlat= 4,urcrnrlat= 50,\
                                                resolution='h', llcrnrlon=45,urcrnrlon=125)
                  map.drawparallels(np.arange(4,50,6),labels=[1,0,0,0],linewidth=0,fontsize=8)
                  map.drawmeridians(np.arange(45,125,8),labels=[0,0,0,1],linewidth=0,fontsize=8)
                  map.drawcoastlines(linewidth=1.5)
                  map.drawcountries(linewidth=1.5)
                  lons, lats = map(lon, lat)
                  map.scatter(lons, lats, marker = 'o', s = 15, color='r')
                  return ax

fig=plt.figure() #figsize=(a,b) may be passed into this
ax1 = plt.subplot(221)
ax2 = plt.subplot(222)
ax3 = plt.subplot(223)
axlist = [ax1,ax2,ax3]
for ax in axlist:
                plot_background(ax)

b1=C.variables['lon'][:]
c1=C.variables['lat'][:]
x,y=np.meshgrid(b1,c1)
CI=[-80,-72,-64,-56,-48,-40,-32,-24,-16,-8,0,8,16,24,32,40,48,56,64]

A = axlist[0].contourf(x,y,ANA1,CI,cmap='gist_ncar',extend='both',zorder=1)
A1 = axlist[0].contour(x,y,ANA1,CI,colors='black',linewidths=1,zorder=2)
axlist[0].clabel(A1, inline=1, fontsize=8,fmt="%i")
axlist[0].set_title('16-10-1977',fontsize=10,pad=10)

B = axlist[1].contourf(x,y,ANA2,CI,cmap='gist_ncar',extend='both',zorder=1)
B1 = axlist[1].contour(x,y,ANA2,CI,colors='black',linewidths=1,zorder=2)
axlist[1].clabel(B1, inline=1, fontsize=8,fmt="%i")
axlist[1].set_title('17-10-1977',fontsize=10,pad=10)

D = axlist[2].contourf(x,y,ANA3,CI,cmap='gist_ncar',extend='both',zorder=1)
D1 = axlist[2].contour(x,y,ANA3,CI,colors='black',linewidths=1,zorder=2)
axlist[2].clabel(D1, inline=1, fontsize=8,fmt="%i")
axlist[2].set_title('18-10-1977',fontsize=10,pad=10)

Z=plt.colorbar(A, ax=[ax1,ax2,ax3],shrink=0.8, aspect=20, fraction=.15,pad=.05)
Z.set_ticks([-80,-72,-64,-56,-48,-40,-32,-24,-16,-8,0,8,16,24,32,40,48,56,64])
Z.set_label('m', labelpad=-40, y=1.05, rotation=0)
plt.suptitle('Geopotential height anomalies for given dates at 850 hPa', fontsize=15) 

Following image is generated after the code runs

Answer 2

这里有一个参考可能很有用，尽管它使用了xarray和cartopy，而不是Matplotlib中的basemap

https://unidata.github.io/MetPy/latest/examples/Four_Panel_Map.html

Xarray在处理多维数据时非常有用。您可以使用sel和slice函数轻松选择给定时间范围内的值

您可以在此处了解更多信息： http://xarray.pydata.org/en/stable/time-series.html