[GRASS-user] Calculate average of 3 highest values for a set of raster maps

[Pietro]

Hi Milton,

On Friday 15 Mar 2013 08:13:16 Milton Cezar Ribeiro wrote:
> Dear all,
>
> Thanks for all replies. I will check it out thoughout the weekend.
> Just clarifying something, I need to do the estimations focusing on
> each pixels, and not for the full map (I was not that clear before).

In GRASS7 you can use the RasterRow or RasterRowIO class of pygrass, something
like:

{{{
#!python
import numpy as np

from grass.pygrass.raster import RasterRow
from grass.pygrass.gis.region import Region
from grass.pygrass.functions import coor2pixel

# define the list of points that you want
POINTS = [(633042.213884, 226283.771107),
          (641197.93621, 225675.891182),
          (644034.709193, 220711.538462),
          (630028.142589, 224764.071295),
          (630002.814259, 215037.992495), ]

NEIGHS = np.array([(-1, 1),
                   (0, 1),
                   (1, 1),
                   (1, 0),
                   (1, -1),
                   (0, -1),
                   (-1, -1),
                   (-1, 0),
                   (-1, 1)])


def get_neighs(point, region, neighs=NEIGHS):
    res = point - neighs
    valid = ((0 <= res.T[0]) * (res.T[0] < region.rows) *
             (0 <= res.T[1]) * (res.T[1] < region.cols))
    return res[valid]


def do_something(rastname, points, mapset=''):
    # get the current region
    region = Region()
    raster1 = RasterRow(rastname, mapset)
    raster1.open('r')
    # instantiate more raster object here if you need
    average = []
    for pnt in points:
        coor = coor2pixel(pnt, region)
        neighs = get_neighs(coor, region)
        sum_ = 0
        for nx, ny in neighs:
            # do something with the neighbour value
            sum_ += raster1[int(nx)][int(ny)]
        average.append(sum_ / len(neighs))
    raster1.close()
    # close the other rasters
    return average

print do_something('elevation', POINTS, mapset='user1')
}}}

it's working on North Carolina mapset... maybe I should add get_neighs into
the pygrass.functions module...
whit the same logic you can open more raster at the same time

> Also as my raster maps are very large (~50,000x60,000 pixels) maybe R
> (at least on w7) maybe will not handle the seven rasters.

I've worked with a 10**5x10**5 without problems [0].

Pietro

[0] http://www.mdpi.com/2220-9964/2/1/201