Thanks Pietro and Paulo I will try all your suggestions. cheers
milton 2013/3/15, Pietro <[email protected]>: > 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 - [email protected] 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 _______________________________________________ grass-user mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/grass-user
