Extracting a variable and averaging through all time slices

[Jacob Alberto Garcia]

Hi!

I have this script that was made from various sources available for wrf-python. What I want to do is to get all time slices and average it through all times, and then plot the final single slice (Averaged). This is for the output from a WRF/Chem run. Any ideas on how I could edit this script to do it? Thanks for any help! 

from netCDF4 import Dataset
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
import cartopy.crs as crs
from cartopy.io.shapereader import Reader
from cartopy.feature import ShapelyFeature
from cartopy.feature import NaturalEarthFeature

from wrf import (to_np, getvar, smooth2d, get_cartopy, cartopy_xlim,
                 cartopy_ylim, latlon_coords)

# Open the NetCDF file
wrf_file = Dataset('/Users/jacob/Desktop/wrfpy/wrfout_d01.nc')

# Get the sea level pressure
pm25 = getvar(wrf_file, "PM2_5_DRY", timeidx=23)[0, :]

# Smooth the sea level pressure since it tends to be noisy near the
# mountains


# Get the latitude and longitude points
lats, lons = latlon_coords(pm25)

# Get the cartopy mapping object
cart_proj = get_cartopy(pm25)

# Create a figure
fig = plt.figure(figsize=(8,6))
# Set the GeoAxes to the projection used by WRF
ax = plt.axes(projection=cart_proj)

# Download and add the states and coastlines



fname = '/Users/jacob/Desktop/wrfpy/phil.shp'

shape_feature = ShapelyFeature(Reader(fname).geometries(),
                                crs.PlateCarree(), facecolor='none')
ax.add_feature(shape_feature, linewidth=2, edgecolor="black")



# Make the contour outlines and filled contours for the smoothed sea level
# pressure.
plt.contour(to_np(lons), to_np(lats), to_np(pm25), 10, colors="black",
            transform=crs.PlateCarree(), linewidths=0.5, linestyles='dotted')
plt.contourf(to_np(lons), to_np(lats), to_np(pm25), 10,
             transform=crs.PlateCarree(),
             cmap=get_cmap("jet"), alpha = 0.3)




# Add a color bar
plt.colorbar(ax=ax, shrink=.98)

# Set the map bounds


# Add the gridlines
ax.gridlines(color="black", linestyle="dotted")

plt.title("Particulate Matter 2.5 microns in Diameter Concentration")

plt.show()