Generating Multiple basemap plots from NetCDF using loop

Here's a reference that might be useful though it uses xarray and cartopy instead of basemap from Matplotlib.

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

Xarray is very useful when dealing with multidimensional data. You can easily select the values within a given time range using sel and slice functions.

you can learn more about it here: http://xarray.pydata.org/en/stable/time-series.html

Thank you Dani56 for the necessary links, I was able to generate the plots. The code for generating multiple Geo-potential height anomaly plots is given below:

#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