With count data, like population, it only makes sense to aggregate using a sum resampling function (like you did with block_reduce).
Using rasterio you can perform a decimated read (setting out_shape) using resampling=Resampling.sum. I see from the docs that it requires GDAL >= 3.1 though (https://rasterio.readthedocs.io/en/latest/api/rasterio.enums.html#rasterio.enums.Resampling)
It can give unexpected results if the file has precomputed overviews, then no matter which resampling method you specify, rasterio won't be able to ignore them (unless built against GDAL >= 3.3 which is very unlikely at the moment).

cheers,
Loïc
________________________________________
From: main@rasterio.groups.io <main@rasterio.groups.io> on behalf of tinak.contact@gmail.com <tinak.contact@gmail.com>
Sent: 17 March 2021 07:06:04
To: main@rasterio.groups.io
Subject: Re: [rasterio] Rasterio read outshape resample does not give correct array

[Edited Message Follows]
[Reason: added gdal option ]

Thank you very much, I found scikit block_reduce<https://urldefense.com/v3/__https://scikit-image.org/docs/dev/api/skimage.measure.html*block-reduce__;Iw!!DOxrgLBm!SYQfUIkp3oT6I4nkx5rhk0VhOIXLuJc03ZvgBB-xB81wfd0QZgYJdbZuqh9ZDbxz4GMJQqQ$> which did the trick, but I was wondering if there are any gdal or rasterio functions who could
do the same. I found from (here<https://urldefense.com/v3/__https://gis.stackexchange.com/questions/152661/downsampling-geotiff-using-summation-gdal-numpy__;!!DOxrgLBm!SYQfUIkp3oT6I4nkx5rhk0VhOIXLuJc03ZvgBB-xB81wfd0QZgYJdbZuqh9ZDbxz3Jieug0$>) and ( here<https://urldefense.com/v3/__https://gdal.org/programs/gdalwarp.html__;!!DOxrgLBm!SYQfUIkp3oT6I4nkx5rhk0VhOIXLuJc03ZvgBB-xB81wfd0QZgYJdbZuqh9ZDbxzYF_dZ3o$>) that
GDAL 3.1 or later does have -r sum or -r average
as interpolation method.
However, if I apply the gdal translate command I am not getting the correct numbers. Could this be a problem with my gdal?


Solution with scikit block_reduce :<https://urldefense.com/v3/__https://scikit-image.org/docs/dev/api/skimage.measure.html*block-reduce__;Iw!!DOxrgLBm!SYQfUIkp3oT6I4nkx5rhk0VhOIXLuJc03ZvgBB-xB81wfd0QZgYJdbZuqh9ZDbxz4GMJQqQ$>

from skimage.measure import block_reduce
import matplotlib.pyplot as plt
factor = 2
with rasterio.open("mwi_ppp_2020_UTM.tif") as dataset:
data = dataset.read(1)
print('population file', data[np.where(data>0)].sum())
# give zero to no dataval
data=np.where(data<=0,0, data)

# reduce by factor
x=block_reduce(data, block_size=(factor,factor), func=np.sum)

# mask out no data vals
x_masked=np.ma.masked_equal(x,0)
plt.imshow(x_masked)
plt.colorbar()
print('population resampled',x_masked.sum())
[cid:attach_0_166D0C1E101B79CD_3895@groups.io]


Using gdal translate -r average

(decrease by factor 10)
import subprocess
subprocess.call("""gdal_translate -tr 0.0083 0.0083 -r average\
/Users/admin/Downloads/mwi_ppp_2020.tif gdal_test.tif""", shell=True)

Population after 206120.61