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

[Milton Cezar Ribeiro]

Thanks Pietro and Paulo
I will try all your suggestions.

cheers

milton

2013/3/15, Pietro <peter.zamb at gmail.com>:
> 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
>


-- 
Miltinho - mcr at rc.unesp.br
Laboratório de Ecologia Espacial e Conservação - LEEC
Depto de Ecologia - UNESP - Rio Claro
Av. 24A, 1515- Bela Vista
13506-900 Rio Claro, SP, Brasil

Fone: +55 19 3526-9647 (office)  19 3526-9680 (lab)
Cel: 19 9853-3220 / 19 9853-5430

Depto Ecologia http://www.rc.unesp.br/ib/ecologia/

PG ECO & BIODIV http://www.rc.unesp.br/ib/ecologia/posbiodiversidade/index.php

CV http://buscatextual.cnpq.br/buscatextual/visualizacv.do?id=K4792988H6&mostrarNroCitacoesISI=true&mostrarNroCitacoesScopus=true

Google citations http://scholar.google.com/citations?user=OWX_2eAAAAAJ