使用循环从 NetCDF 生成多个底图答案

【问题标题】：Generating Multiple basemap plots from NetCDF using loop使用循环从 NetCDF 生成多个底图
【发布时间】：2020-04-26 11:28:47
【问题描述】：

我需要生成多个位势高度异常图（数据取自 NCEP NCAR 再分析），为此我使用 python netCDF4 和底图包。我编写的代码生成单个图像，但想修改它以生成多个图。生成单张图片的代码如下：

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

我有 2 个问题 1. 如何存储多个 NetCDF 文件中日期 177 到 180 和 2.如何使用存储的数据生成多个底图，即从177到180的所有日期[日期177代表27-06-2010等]

【问题讨论】：

标签： python matplotlib-basemap

【解决方案1】：

这里有一个参考，虽然它使用 xarray 和 cartopy 而不是 Matplotlib 中的底图，但它可能很有用。

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

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

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

【讨论】：

【解决方案2】：

感谢 Dani56 提供必要的链接，我能够生成图。生成多个地势高度异常图的代码如下：

#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

【讨论】：