Apply mask and compute statistics with GDAL

I've constructed a VRT that includes a band of data and a mask (0/1 representing water/land) with:

``````ds = gdal.Translate(masked_vrt, combined_vrt, bandList=, maskBand=2, format='VRT')
``````

Where combined VRT has band 1 of data and band 2 of 0/1.

I would like to use the python GDAL bindings with this VRT to compute the statistics of the data after applying the mask but `ds=gdal.Open('masked.vrt').GetRasterBand(1).GetStatistics(0,1)` gives only the stats of the original, maskless, data (band 1).

I can get the result I'm looking for like this:

``````ds = gdal.Open('masked.vrt')
``````

But reading in the arrays of course loses the efficiency of using the VRT.

Is there some way to use GDAL directly to apply the mask and compute the stats?

Here is a MWE:

``````import numpy as np
from osgeo import gdal, gdal_array
data1 = np.ones([5, 5]) * 2
data2 = np.ones([5, 5]) * 5
data  = np.vstack([data1, data2])

## dummy geotransform
gt = [-78.38, 2.7e-04, 0, 37.33, 0, -2.7e-04]

data_ds = gdal_array.OpenArray(data)
data_ds.SetGeoTransform(gt)

## write datasets
gdal.Translate('data', data_ds, outputType=gdal.GDT_Float64)
gdal.BuildVRT('data.vrt', 'data')