Python bindings: Add Dataset.ReadAsMaskedArray

[dbaston]

Resolves #12005

[rouault]

Perhaps indicate in the doc that when downsampling the mask resampling algorithm is Mode (ie majority) ?

[dbaston]

Documented this and also changed behavior of Band.ReadAsMaskedArray to match.

[dbaston]

Hmm, I'm not sure how my tests were passing locally before. Right now I'm puzzled by

(Pdb) band.GetMaskBand().ReadAsArray()
array([[  0, 255],
       [  0,   0]], dtype=uint8)
(Pdb) band.GetMaskBand().ReadAsArray(buf_xsize=1, buf_ysize=1, resample_alg=gdal.GRIORA_Mode)
array([[255]], dtype=uint8)

[rouault]

Right now I'm puzzled by

I believe this is due to: https://github.com/OSGeo/gdal/blob/master/gcore/overview.cpp#L4453 / https://github.com/OSGeo/gdal/blob/master/gcore/overview.cpp#L5798

[dbaston]

I don't think GRA_Mode downscaling works for NoData mask bands. If I'm following the code right, it performs a downscaled read from the parent band (using GRA_Mode, and excluding NoData pixels) and then sets all valid pixels to 255. Which is different from performing a full resolution read on the parent, setting valid pixels to 255, and then downscaling.

[rouault]

I don't think GRA_Mode downscaling works for NoData mask bands.

yes the logic I pointed out where the mask band of the mask band is the mask band itself should probably be modified for Mode, and just get the GetMaskBand() of the mask band, which would be a AllValid mask band

[dbaston]

I believe this is due to: https://github.com/OSGeo/gdal/blob/master/gcore/overview.cpp#L4453 / https://github.com/OSGeo/gdal/blob/master/gcore/overview.cpp#L5798

I don't see either of these lines hit in this case. The issue I see is

gdal/gcore/gdalnodatamaskband.cpp

Lines 404 to 407 in bfb3825

	eErr = m_poParent->RasterIO(
	GF_Read, nXOff, nYOff, nXSize, nYSize, pTemp, nBufXSize, nBufYSize,
	eWrkDT, nWrkDTSize, static_cast<GSpacing>(nBufXSize) * nWrkDTSize,
	psExtraArg);

where a downscaled read is performed from the parent band, and the result reclassified to 0/255. To get a correct result for GRIORA_Mode, I think we need to read from the parent band at full resolution, reclassify, and then downscale.

All this said, the current behavior of GRIORA_Mode is essentially GRIORA_Max and not necessarily undesirable as a default.

[dbaston]

@rouault any thoughts on me marking test_numpy_rw_masked_array_mask_downsampling as xfail, then proceeding along these lines in a separate PR?

I think we need to read from the parent band at full resolution, reclassify, and then downscale

[rouault]

then proceeding along these lines in a separate PR?

I think we need to read from the parent band at full resolution, reclassify, and then downscale

Only for mode ? Otherwise I'm slightly concerned that if using a different downscaling algorithm for the mask than the actualband, there could be situations where this would lead to the downscaled mask indicating a pixel as valid whereas the actual downscaled band would have determined it not to be valid. In particular, for cubic or lanczos, there are subtelties regarding that determination (cf overview.cpp lines 3447-3471 (horizontal pass), 3654-3684 (vertical pass))

[dbaston]

Only for mode ?

It's hard for me to picture what cubic resampling of a binary mask should yield, so it wouldn't bother me to keep the status quo there.

I'm approaching this from the expectation that resampling a mask band should yield the same results as resampling an ordinary band with the same values (0 or 255). Is that the wrong way to look at it?

[rouault]

Is that the wrong way to look at it?

I'm not sure what the "right" way is. I'm just underlying that the overview resampling code (also used by RasterIO() for non-nearest resampling) uses the mask band as its own mask band, ignoring 0 values. Both approaches (current GDAL one or the one you suggest) have their pros and cons depending use cases (if you want to be "optimistic" or "pessimistic" regarding pixel validity). Perhaps we should have a way for the user to specify which it prefers, but that would further complicates things...

[dbaston]

I've exposed a mask_resample_alg argument in both ReadAsMaskedArray() methods, which I think is a resonable addition: I can imagine a user wanting to read pixels with GRIORA_Average and have the mask resampled with GRIORA_Mode, GRIORA_Min, or GRIORA_Max (I realize these last two don't currently exist).

This also sidesteps questions about the behavior of those algorithms for now.

[coveralls]

coverage: 70.754% (+0.01%) from 70.74%
when pulling 4f495bc on dbaston:python-dataset-read-as-masked-array
into 4aa3135 on OSGeo:master.