There are a bunch of open source tools beyond the GDAL NetCDF driver you could use. These tools are specialized for your outlined problem area of weather/climate statistics.
- You could have a look "Climate Data Operators" to work with NetCDf /GRIB encoded data of the ERA-5 ECMWF reanalysis data. CDO is focused on "climate data" with as well defined interface. This interface addresses an international standard with a specialized NetCDF/GRIB data configuration for provenience, dimensions (spatial, time) and measured/calculated records (pressure, temperature ...). AFAIK most ECMWF data products follow this standard due to the need of international exchange of data sets in the context of (European Center For) Medium-Range Weather Forcast.
The toolkit is very fast and has an operator set for statistics -- section 1.6.4 -- page 89 ff.. For every latitude is an
operator
over all longitudes defined. Operators for zonal statistics are:
cdo [zonemin,zonmax,zonmean,zonmean,zonvar,zonstd] input-file output-file
cdo zonepctl,num input-file output-file
- Another toolkit is "NetCDF operators". To get data out of the NetCDF result set, you could have a look at "How can I extract data from a NetCDF file for a specific location?". NCO has also a set of statistic operators -- section 4.6 -- page 254 ff.. The syntax for an aggregation (min, max, avg, stddev, etc.) of data records over an dimension (spatial i.e. longitude, latitude and time...) is more complex, because you can use a free NetCDF format to store the matter.
nces [-3] [-4] [-5] [-6] [-7] [-A] [-C] [-c] [--cb y1,y2,m1,m2,tpd]
[--cnk_byt sz_byt] [--cnk_csh sz_byt] [--cnk_dmn nm,sz_lmn]
[--cnk_map map] [--cnk_min sz_byt] [--cnk_plc plc] [--cnk_scl sz_lmn]
[-D dbg] [-d dim,[min][,[max][,[stride]]] [-F]
[-G gpe_dsc] [-g grp[,...]] [--glb ...]
[-h] [--hdf] [--hdr_pad nbr] [--hpss]
[-L dfl_lvl] [-l path] [-n loop]
[--no_cll_msr] [--no_frm_trm] [--no_tmp_fl] [--nsm_fl|grp] [--nsm_sfx sfx]
[-O] [-o output-
file] [-p path] [--ppc ...] [-R] [-r] [--ram_all] [--rth_dbl|flt]
[-t thr_nbr] [--unn] [-v var[,...]] [-w wgt] [-X ...] [-x] [-y op_typ]
[input-files] [output-file]
- For completeness I will mention the NetCDF Command Language. You could write scripts in NCL and have predefined functions for statistics too. NCL produces terrific predefined maps. These maps are good, if you need an orientation and the use cases are well documented. The example page for zonal statistics provides eight use cases. NCL example for a zonal average temperature (www.ncl.ucar.edu).

;*************************************************
; zonal_3.ncl
;
; Concepts illustrated:
; - Attaching a zonal means plot to a cylindrical equidistant contour plot
;************************************************
;
; These files are loaded by default in NCL V6.2.0 and newer
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
;************************************************
begin
;************************************************
; variable and file handling
;************************************************
in = addfile("83.nc","r")
ts = in->TS ; select variable to ave
;************************************************
; plotting
;************************************************
wks = gsn_open_wks("png","zonal") ; send graphics to PNG file
res = True ; make plot mods
res@tiMainString = "Zonal Average" ; Title for the plot
res@gsnZonalMean = True ; put zonal on side
res@gsnZonalMeanXMinF = 240. ; set minimum X-axis value for zonal mean plot
res@gsnZonalMeanXMaxF = 315. ; set maximum X-axis value for zonal mean plot
res@gsnZonalMeanYRefLine = 273.15 ; set reference line X-axis value
plot=gsn_csm_contour_map(wks,ts(0,:,:),res) ; plot temp contours
;***********************************************
end