Re: Rasterio read outshape resample does not give correct array


tinak.contact@...
 
Edited

Thank you very much, I found  scikit block_reduce 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) and ( here) 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 :

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())



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

Join main@rasterio.groups.io to automatically receive all group messages.