[GRASS-dev] Extracting coordinates of RasterNumpy objects (pygrass)

[Pietro]

Hi Pierre,

On Tue, Dec 18, 2012 at 1:50 AM, Pierre Roudier
<pierre.roudier at gmail.com> wrote:>
> - Is there a method to create a set of points (class Point) and write
> it into the GRASS database?

yes, but at the moment supports only 2D features... and I don't
understand why... :-)

it was supposed to work with:

# python source code
import numpy as np
from grass.pygrass.raster import RasterNumpy
from grass.pygrass.vector import VectorTopo
from grass.pygrass.vector.geometry import Point
from grass.pygrass.region import Region
from grass.pygrass.functions import pixel2coor

elev = RasterNumpy("elevation")
elev.open()
a = elev > 144

reg = Region()

nonzero = a.nonzero()
coords = np.array([pixel2coor(pixel, reg) for pixel in
np.column_stack(nonzero)])
points = np.rec.fromarrays([coords.T[0], coords.T[1], elev[nonzero]])


new = VectorTopo('newvect')
new.open("w", overwrite=True)
for pnt in points:
    new.write(Point(*pnt))
new.close()
new.build()
# end source code

actually, even if the point is a 3D point, the write method of the
Vector class it is writing only a 2D point... I miss something but I
don't know where...
I hope to fix this as soon as possible...

> - I would like to use this {x,y,z} set of points in v.surf.bspline.
> don't suppose I got another choice but write my numpy array into the
> GRASS database first right?

I don't think is in the database... I think it's enough if you write
the points in the GRASS topology file... And of course I think you
must do it, if you want to use the v.surf.bspline module...

Pietro