How to resize a raster SUMMING the values of pixels in QGIS?

Question

I would like to be able to resize a raster such that the values of the resized pixels are based on the SUMS of the original raster pixels, rather than the averages. The output should be generated thus: each new pixel value equals the sum of any old pixel values contained fully within its "zone" of coverage, plus the corresponding fraction of the values of any pixels which are partially within its "zone" of coverage.

This would be very useful for resizing rasters containing certain types of data for which averaging the values is not necessarily appropriate, and preserving exactly the same total sum of pixel values may be important--population counts, for example. Is there a straightforward way to do this in QGIS?

Answer 1

After further testing, I've discovered that r.resamp.stats actually does NOT perform this operation correctly, at least for sums. For whatever reason, I have been unable to perform the summing rescale without a "loss" of mass between 25% and 1400%.

I've coded a solution in Python, which for the target resolution I'm working with (1 degree) executes quickly enough. As long as the source array is in double-precision float format, the "loss" is effectively zero (1e-15 % or so).

The code below can be applied to a raster which has already been loaded as a numpy array of any data type.

import numpy as np

def summingit(arrayin,areas,rowsin,colsin,rowsout,colsout,curtargr,curtargc):
     firstinr = float(curtargr)*(float(rowsin)/float(rowsout))
     firstinrw = 1 - firstinr + int(firstinr)
     firstinr = int(firstinr)

     lastinr = (float(curtargr) + 1.)*(float(rowsin)/float(rowsout))
     lastinrw = lastinr - int(lastinr)
     lastinr = int(lastinr)

     firstinc = float(curtargc)*float((colsin)/float(colsout))
     firstincw = 1 - firstinc + int(firstinc)
     firstinc = int(firstinc)

     lastinc = (float(curtargc) + 1.)*(float(colsin)/float(colsout))
     lastincw = lastinc - int(lastinc)
     lastinc = int(lastinc)

     rowslength = int(lastinr - firstinr)
     if lastinrw > 0.:
         rowslength += 1
         if rowslength == 1:
             rowweights = [lastinrw + firstinrw - 1.]
         else:
             rowweights = [firstinrw]
             rowweights.extend(np.ones(rowslength - 2).tolist())
             rowweights.extend([lastinrw])
     else:
         lastinr -= 1
         if rowslength == 1:
             rowweights = [firstinrw]
         else:
             rowweights = [firstinrw]
             rowweights.extend(np.ones(rowslength - 1).tolist())

     rowinds = np.linspace(firstinr,lastinr,rowslength,dtype = 'uint32')

     colslength = int(lastinc - firstinc)
     if lastincw > 0.:
         colslength += 1
         if colslength == 1:
             colweights = [lastincw + firstincw - 1.]
         else:
             colweights = [firstincw]
             colweights.extend(np.ones(colslength - 2).tolist())
             colweights.extend([lastincw])
     else:
         lastinc -= 1
         if colslength == 1:
             colweights = [firstincw]
         else:
             colweights = [firstincw]
             colweights.extend(np.ones(colslength - 1).tolist())

     colinds = np.linspace(firstinc,lastinc,colslength,dtype = 'uint32')

     rowinds = np.array(([rowinds],)*colslength).transpose().flatten()
     colinds = np.array(([colinds],)*rowslength).flatten()
     weights = np.outer(rowweights,colweights).flatten()

     return np.sum(arrayin[rowinds,colinds]*areaweights*weights)

curarray = np.array(curarray,dtype='float64')
rowinds = np.array([range(0,180),]*360).transpose().flatten()
colinds = np.array([range(0,360),]*180).flatten()
newarraytest = map(lambda rows,cols: summingit(curarray,2160,4320,180,360,rows,cols),rowinds,colinds)
newarraytest = np.reshape(newarraytest,(180,360))

In current form the routine will only work for scaling DOWN (i.e., the same operation that r.resamp.stats is supposed to do), though the base function "summingit" could easily be rolled into a different routine that scales UP.

Hope this may be helpful to someone in the future!

Answer 2

Never done it before, but I think that the GRASS algorithm r.resamp.stats should be able to do what you're looking for. You can run the algorithm through QGIS's processing toolbox. The documentation for r.resamp.stats is here: https://grass.osgeo.org/grass64/manuals/r.resamp.stats.html