使用循环从NetCDF生成多幅基图图谱

问题描述 投票：0回答：1

我需要生成多张地势高度异常图（数据来自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.如何使用存储的数据生成多个basemap图，即177到180的所有日期[177日期代表27-06-2010，以此类推]。

python matplotlib-basemap
1个回答
0
投票

这里有一个参考文献，可能会有用，尽管它使用xarray和cartopy而不是Matplotlib的basemap。

https:/unidata.github.ioMetPylatestexamplesFour_Panel_Map.html。

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

你可以在这里了解更多关于它的信息。http:/xarray.pydata.orgenstabletime-series.html。

0
投票

谢谢你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)

以下是代码运行后生成的图片

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.