It was spotted that when summing a float32 layer, the answer didn't match if you summed the entire layer directly using numpy. Worse, the result varied slightly depending on the size of YSTEP.
The root cause was that the accumulator in sum was a float64, so both types ended up being involved in the equation, so there was mixing of precisions that lead to different results.
Now if you sum a float32 layer we force maximum precision by moving everything to float64. This makes sense as the final result was always float64 (a debatable position, perhaps it should match the layer type, but that's for a discussion for another podcast). This also matches the result in QGIS, which is probably useful also.
It was spotted that when summing a float32 layer, the answer didn't match if you summed the entire layer directly using numpy. Worse, the result varied slightly depending on the size of YSTEP.
The root cause was that the accumulator in sum was a float64, so both types ended up being involved in the equation, so there was mixing of precisions that lead to different results.
Now if you sum a float32 layer we force maximum precision by moving everything to float64. This makes sense as the final result was always float64 (a debatable position, perhaps it should match the layer type, but that's for a discussion for another podcast). This also matches the result in QGIS, which is probably useful also.