Генерация нескольких графиков базовой карты из NetCDF с использованием цикла - PullRequest
Новая
       10

Генерация нескольких графиков базовой карты из NetCDF с использованием цикла

0 голосов
/ 26 апреля 2020

Мне нужно сгенерировать несколько графиков геопотенциальных аномалий высоты (данные взяты из реанализа 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 и т. д.]

1 Ответ

0 голосов
/ 07 мая 2020

Вот ссылка, которая может быть полезна, хотя она использует 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

...