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=[1], 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')
data = ds.GetRasterBand(1).ReadAsArray()
mask = ds.GetRasterBand(1).GetMaskBand().ReadAsArray()
masked_data = data[mask>0].mean()

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
mask1 = np.zeros([5, 5])
mask2 = np.ones([5, 5])
data  = np.vstack([data1, data2])
mask  = np.vstack([mask1, mask2])

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

data_ds = gdal_array.OpenArray(data)

mask_ds = gdal_array.OpenArray(mask)

## write datasets
gdal.Translate('data', data_ds, outputType=gdal.GDT_Float64)
gdal.BuildVRT('data.vrt', 'data')
gdal.Translate('mask', mask_ds, outputType=gdal.GDT_Float64)
gdal.BuildVRT('mask.vrt', 'mask')
gdal.BuildVRT('combined.vrt', ['data.vrt', 'mask.vrt'], separate=True)

ds = gdal.Translate('applied_mask.vrt', 'combined.vrt', bandList=[1], maskBand=2, format='VRT')
ds1 = gdal.Open('data.vrt', gdal.GA_ReadOnly)
print ('Original Statistics:', ds1.GetRasterBand(1).GetStatistics(0, 1))
print ('Masked Statistics:', ds.GetRasterBand(1).GetStatistics(0,1))
ds = None; ds1 = None

If working, original and masked would differ, of course.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.