I have a Land Cover classification map that I created from 2m resolution multispectral imagery and I am trying to downsample the results to a 20m resolution class map.

The climate in question is arid so there are many areas with sparse shrubs or trees. As such, one 2m resolution class may have a couple of pixels for a single shrub. I am not interested in where each individual tree or bush is, but rather where there are areas of high concentration of a given class (i.e. I want to create a 20m resolution class map that will aggregate the high resolution classes into classes such as "sparse shrubs" or "dense shrubs" depending on the density of 2m x 2m shrub classes found in the 20m x 20m grid square)

Using r.mapcalc in GRASS GIS I know you can refer to neighboring cells using the format map[1,-2] and there are many useful functions available in r.mapcalc. However, my problem is that when I am downsampling to ~20m resolution from 2m resolution there are ~100 neighbors to analyze and I would have to address each one specifically as there is no way that I have found to nest a for loop within a call to r.mapcalc.

Does anyone have suggestions on a way to gather statistics of the cells surrounding a given map cell, and change the cell in question based on its neighbors in the way I describe?

  • If you think that is necessary to use a loop with r.mapcal then you need some script language. In Linux, you can use bash script and, in Windows, a bat file. Another option is using grass.script with python where the sintaxis of r.mapcalc changes to 'grass.mapcalc'. However, I would use GDAL-Python. You can see this approach: gis.stackexchange.com/questions/72150/…. The 3x3 block can be changed for a 10x10 block (100 neighbors).
    – xunilk
    Commented May 5, 2015 at 18:15
  • @xunilk: I saw some examples of people using scripts to run r.mapcalc multiple times, but it didn't seem like that afforded the flexibility I needed when it came to what r.mapcalc was doing on a specific pass. However, your suggestion of using GDAL has been very helpful so far! Thank you
    – G0neSailin
    Commented May 6, 2015 at 19:19
  • Don't you want to be using r.resample or r.neighbors with method=mode?
    – Micha
    Commented May 7, 2015 at 5:33
  • @Micha: That doesn't accomplish my purpose and I do not just want to know the most common value within my window. Instead, I would like to be able to tabulate how many of each class type are within the window. For example, 80% of the cells are sand and 20% of the cells have a class type of shrub. Then in my aggregrate 20m x 20m class map I can equate that 20m resolution cell to "shrub, sparse". Your suggestion would only yield "sand" which is of less interest to me.
    – G0neSailin
    Commented May 7, 2015 at 17:01
  • I guess I didn't quite understand your question. If you want to do any kind of weighted average of the Landcover classifications then you'll loose the classifications. Do you already have in mind a set of new classifications such as "sparse vegetation", "dense vegetation" with the values that you expect, as percents of the original classes?
    – Micha
    Commented May 8, 2015 at 9:37

1 Answer 1


As @xunilk suggested, I found GDAL to be the tool I was looking for. I also used Counter to create a form of histogram that helped me complete my analysis.

for i in xrange(0,iRange,resampleInterval):
    for j in xrange(0,jRange,resampleInterval):
        scanline = band.ReadRaster(i, j, windowSize, windowSize, windowSize, windowSize, band.DataType)
        neighbors = struct.unpack(fmtTypes[BandType] * windowSizeSquared, scanline)
        neighborData = Counter(neighbors)
        statsTuple = neighborData.most_common() # Returns all unique items and their counts [(4,101),(1,50)]
        numNeighbors = sum(neighborData.values()) # sum(neighborData.values()) yields number of neighbors
        # findClass assigns classID based on percent of total area each class makes up
        classID = findClass(statsTuple, numNeighbors, classID, 'Roads')
        classID = findClass(statsTuple, numNeighbors, classID, 'Water')

        classID = findClass(statsTuple, numNeighbors, classID, 'Forest')
        classID = findClass(statsTuple, numNeighbors, classID, 'Shrub')
        # If mixture does not fit into one of the above buckets, set class equal to the mode of the window
        if not classID:
            classID = statsTuple[0][0] # Set classID equal to mode of neighbors
        classArray.append((i+1,j+1,classID)) #save x,y coordinate with classID

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.