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