<p>谢谢你提供必要的链接,我能够生成情节。生成多个地质位高度异常图的代码如下所示:</p>
<pre><code>#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)
</code></pre>
<p><a href="https://i.stack.imgur.com/LEG0Q.png" rel="nofollow noreferrer">Following image is generated after the code runs</a></p>