Re: Rasterio read outshape resample does not give correct array
Loïc Dutrieux
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
|
|