0

I have a population grid and want to decrease the resolution by a factor 2. I would like to get a raster that has 2 times less rows and columns, a resolution that is 2 times as large as the original one, while retaining the same number of population.

This should be a straightforward operation, but the population numbers change incorrectly if I follow the documentation. When the resampled array is written to file, the numbers are incorrect.

import rasterio 
import numpy as np 
import matplotlib.pyplot as plt
from rasterio.enums import Resampling
 
 
factor = 2

with rasterio.open("mwi_ppp_2020.tif") as dataset:
    
    arr = dataset.read(1,masked=True)
    #Show total population of Malawi 
    print(arr.sum())

    data = dataset.read(1, masked=True,resampling=Resampling.bilinear,
        out_shape=(
            int(dataset.height * factor),
            int(dataset.width * factor)))

     #Show total population of array 
    print(data.sum())

     # scale image transform
    transform = dataset.transform * dataset.transform.scale(
            (dataset.width / data.shape[1]),
            (dataset.height / data.shape[0])
                                )
    

enter image description here

Shouldn't the resample method yield an array with resampled cell values to get the same number of total population?

5
  • I don't think so? Looks like you're increasing the total number of pixels then you're summing all of their values. What are you trying to do exactly?
    – mikewatt
    Mar 16, 2021 at 19:43
  • Thanks mikewatt. I would like to get raster that has 2 times less rows and columns, a resolution that is 2 times as large as the original one, while retaining the same number of population
    – Tina K
    Mar 16, 2021 at 20:11
  • I was asking more why you want lower resolution data instead of using the raster as-is. But using "sum" as the resampling method might get you close: rasterio.readthedocs.io/en/latest/api/…
    – mikewatt
    Mar 16, 2021 at 21:17
  • The numbers will be correct for the specified resampling algorithm (bilinear interpolation), however they're not what you are wanting. I suspect you want to sum which is only available if your rasterio is using a very recent GDAL (>= 3.1). An alternative is the scikit block_reduce function.
    – user2856
    Mar 17, 2021 at 1:00
  • thank you so much for your replies! I solved the problem by both using a newer version of gdal and alternatively using scikit block_reduce
    – Tina K
    Mar 18, 2021 at 5:26

3 Answers 3

0

I solved the issue using help in the comments. I needed gdal>3.1 for the sum aggregation/interpolation method (bilinear is wrong in this case).

The gdal code is

gdalwarp -tr 0.0083 0.0083 -r sum  -srcnodata -9999 -dstnodata -9999 mwi_ppp_2020.tif  outtest.tif

scikit block_reduce was also working as alternative (with factor 10 in this example), but you have to very careful with the nodata entries and the affine transform. So gdalwarp is the better alternative.

enter image description here

0

If gdal is used, then it's also possible to use the gdal python bindings.

from osgeo import gdal

# Input and output file paths
input_file = "mwi_ppp_2020.tif"
output_file = "outtest.tif"

# Resampling method
resample_alg = "sum"

# Warp options
options = gdal.WarpOptions(
    format="GTiff",
    xRes=0.0083,
    yRes=0.0083,
    srcNodata=-9999,
    dstNodata=-9999,
    resampleAlg=resample_alg
    )

# Perform the warp operation
gdal.Warp(output_file, input_file, options=options)

gdal.Warp() and gdal.WarpOptions() docs.

-1

Try:

int(dataset.height / factor),
int(dataset.width / factor)))
1
  • Sadly, this does not work. It reduces the population count to about 4million and does not calculate the correct cell values
    – Tina K
    Mar 16, 2021 at 20:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.