3

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?

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

  • From the documentation, this algorithm does appear to do what I need. Would be nice if it would generalize to also increase resolution according to the same logic, but for the moment I only need to decrease resolution, so should do the trick. I'll try it and report back on how it works. – Matt D. Jul 19 '16 at 9:11
  • great! glad to hear it – jbukoski Jul 20 '16 at 15:04
  • Sorry to de-accept your answer, but after further investigation it turns out that r.resamp.stats doesn't work like we both thought it should for sums. I've posted the python code for my solution. – Matt D. Aug 1 '16 at 18:48
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!

  • Could you please provide an example of how a raster is imported in the required format to run your code? – Jackk Mar 29 at 10:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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