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

Pierre_Roudier · December 17, 2012, 6:48am

Hi,

I am using (with enthusiasm!) Pietro’s pygrass library to develop a raster module. I am using Numpy/Scipy as my working horse, so I manipulate a lot of the RasterNumpy objects that have been introduced with pygrass.

In a specific step, I am identifying pixels using a test (this would be similar to my_array < 100). I would like to extract the points that satisfy the test, and access not only their values and index in the array, but I would also go back to their coordinates to extract them as (x, y, z).

Is there a way to do this using RasterNumpy and pygrass?

Cheers,

Pierre

Pietro2 · December 17, 2012, 8:45am

Hi Pierre,

On Mon, Dec 17, 2012 at 7:48 AM, Pierre Roudier
<pierre.roudier@gmail.com> wrote:

Hi,

I am using (with enthusiasm!) Pietro's pygrass library to develop a raster
module. I am using Numpy/Scipy as my working horse, so I manipulate a lot of
the RasterNumpy objects that have been introduced with pygrass.

I'm really happy that you are using the pygrass library!

In a specific step, I am identifying pixels using a test (this would be
similar to my_array < 100). I would like to extract the points that satisfy
the test, and access not only their values and index in the array, but I
would also go back to their coordinates to extract them as (x, y, z).

Is there a way to do this using RasterNumpy and pygrass?

yes, you should use the numpy stuff, something like:

import numpy as np
a = np.arange(15).reshape(3, 5)
a

array([[ 0, 1, 2, 3, 4],
[ 5, 6, 7, 8, 9],
[10, 11, 12, 13, 14]])

b = a>7
b1

array([[False, False, False, False, False],
[False, False, False, True, True],
[ True, True, True, True, True]], dtype=bool)

b.nonzero() # return two array with x and y and z if the array is 3D

(array([1, 1, 2, 2, 2, 2, 2]), array([3, 4, 0, 1, 2, 3, 4]))

Please let me know if you find something that is not clear, or is not
working properly...

All the best!

Pietro

Pierre_Roudier · December 17, 2012, 8:44pm

Thanks Pietro,

Yes that answers the question, but just partly: I was actually
wondering whether it would be possible to get the extracted x and y
coordinates in geographic space (as opposed to the Numpy array space)?

Thanks again for the great work on pygrass,

P

2012/12/17 Pietro <peter.zamb@gmail.com>:

Hi Pierre,

On Mon, Dec 17, 2012 at 7:48 AM, Pierre Roudier
<pierre.roudier@gmail.com> wrote:

Hi,

I am using (with enthusiasm!) Pietro's pygrass library to develop a raster
module. I am using Numpy/Scipy as my working horse, so I manipulate a lot of
the RasterNumpy objects that have been introduced with pygrass.

I'm really happy that you are using the pygrass library!

In a specific step, I am identifying pixels using a test (this would be
similar to my_array < 100). I would like to extract the points that satisfy
the test, and access not only their values and index in the array, but I
would also go back to their coordinates to extract them as (x, y, z).

Is there a way to do this using RasterNumpy and pygrass?

yes, you should use the numpy stuff, something like:

import numpy as np
a = np.arange(15).reshape(3, 5)
a

array([[ 0, 1, 2, 3, 4],
       [ 5, 6, 7, 8, 9],
       [10, 11, 12, 13, 14]])

b = a>7
b1

array([[False, False, False, False, False],
       [False, False, False, True, True],
       [ True, True, True, True, True]], dtype=bool)

b.nonzero() # return two array with x and y and z if the array is 3D

(array([1, 1, 2, 2, 2, 2, 2]), array([3, 4, 0, 1, 2, 3, 4]))

Please let me know if you find something that is not clear, or is not
working properly...

All the best!

Pietro

--
Scientist
Landcare Research, New Zealand

Soeren_Gebbert2 · December 17, 2012, 9:11pm

Hi Pierre,
in pygrass you can find the function pixel2coor((col, row), region)
that will transform your column and row into a (northing, easting)
tuple. It is defined in pygrass/functions.py . Internally
Rast_row_to_northing() and Rast_col_to_easting() from the GRASS raster
library are used.

Best regards
Soeren

2012/12/17 Pierre Roudier <pierre.roudier@gmail.com>:

Hi,

I am using (with enthusiasm!) Pietro's pygrass library to develop a raster
module. I am using Numpy/Scipy as my working horse, so I manipulate a lot of
the RasterNumpy objects that have been introduced with pygrass.

In a specific step, I am identifying pixels using a test (this would be
similar to my_array < 100). I would like to extract the points that satisfy
the test, and access not only their values and index in the array, but I
would also go back to their coordinates to extract them as (x, y, z).

Is there a way to do this using RasterNumpy and pygrass?

Cheers,

Pierre

_______________________________________________
grass-dev mailing list
grass-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev

Pietro2 · December 17, 2012, 9:24pm

Hi Pierre,

On Mon, Dec 17, 2012 at 9:44 PM, Pierre Roudier
<pierre.roudier@gmail.com> wrote:

Thanks Pietro,

Yes that answers the question, but just partly: I was actually
wondering whether it would be possible to get the extracted x and y
coordinates in geographic space (as opposed to the Numpy array space)?

aah sorry, I forgot to say that you can use pixel2coord... see the
example below:

from grass import pygrass
from pygrass.raster import RasterNumpy
from pygrass.region import Region
from pygrass.functions import pixel2coor
elev = RasterNumpy("elevation")
elev.open()
a = elev > 144
reg = Region()
xcoords, ycoords = a.nonzero()
for pixel in zip(list(xcoords), list(ycoords))[:5]:

... print pixel
...
(0, 6)
(0, 7)
(0, 8)
(0, 9)
(0, 10)

for pixel in zip(list(xcoords), list(ycoords))[:5]:

... print pixel2coor(pixel, reg)
...
(228440.0, 630000.0)
(228430.0, 630000.0)
(228420.0, 630000.0)
(228410.0, 630000.0)
(228400.0, 630000.0)

Best regards

Pietro

Pierre_Roudier · December 17, 2012, 9:41pm

Thanks Pietro and Soeren,

That's really great!

I am looking forward to use this with v.surf.bspline to do on the fly
interpolation (interpolation with SciPy is quite slow and memory
hungry). I am coding an iterative process that calls the interpolation
routine a lot of times (~10 to 40 interpolation steps), so I would
like to make the most of GRASS interpolation power.

Cheers,

Pierre

2012/12/18 Pietro <peter.zamb@gmail.com>:

Hi Pierre,

On Mon, Dec 17, 2012 at 9:44 PM, Pierre Roudier
<pierre.roudier@gmail.com> wrote:

Thanks Pietro,

Yes that answers the question, but just partly: I was actually
wondering whether it would be possible to get the extracted x and y
coordinates in geographic space (as opposed to the Numpy array space)?

aah sorry, I forgot to say that you can use pixel2coord... see the
example below:

from grass import pygrass
from pygrass.raster import RasterNumpy
from pygrass.region import Region
from pygrass.functions import pixel2coor
elev = RasterNumpy("elevation")
elev.open()
a = elev > 144
reg = Region()
xcoords, ycoords = a.nonzero()
for pixel in zip(list(xcoords), list(ycoords))[:5]:

... print pixel
...
(0, 6)
(0, 7)
(0, 8)
(0, 9)
(0, 10)

for pixel in zip(list(xcoords), list(ycoords))[:5]:

... print pixel2coor(pixel, reg)
...
(228440.0, 630000.0)
(228430.0, 630000.0)
(228420.0, 630000.0)
(228410.0, 630000.0)
(228400.0, 630000.0)

Best regards

Pietro

--
Scientist
Landcare Research, New Zealand

Pierre_Roudier · December 18, 2012, 12:50am

Sorry Pietro - more questions as I'm working my way through this new tool!

- As a follow-up from my last question, I am trying to convert the
pixels I have identified in my raster to a vector layer. While I
manage to create a Point() from scratch, I can't find out how I could
append coordinates to add more points to it. Using the previous
example:

from grass import pygrass
from pygrass.raster import RasterNumpy
from pygrass.region import Region
from pygrass.functions import pixel2coor
elev = RasterNumpy("elevation")
elev.open()
a = elev > 144
reg = Region()
xcoords, ycoords = a.nonzero()
# Switch to geographic space
coords = [pixel2coor(pixel, reg) for pixel in zip(list(xcoords), list(ycoords))]
# Add attribute column
xyz = np.column_stack((np.asarray(coords), z))
from pygrass.vector import *
p = pygrass.vector.geometry.Point() # Works
p = pygrass.vector.geometry.Point(x = xyz[:,0], y = xyz[:,1], z=xyz[:,2]) # Fails

- Is there a method to create a set of points (class Point) and write
it into the GRASS database?

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

Thanks again for your work,

Pierre

2012/12/18 Pierre Roudier <pierre.roudier@gmail.com>:

Thanks Pietro and Soeren,

That's really great!

I am looking forward to use this with v.surf.bspline to do on the fly
interpolation (interpolation with SciPy is quite slow and memory
hungry). I am coding an iterative process that calls the interpolation
routine a lot of times (~10 to 40 interpolation steps), so I would
like to make the most of GRASS interpolation power.

Cheers,

Pierre

2012/12/18 Pietro <peter.zamb@gmail.com>:

Hi Pierre,

On Mon, Dec 17, 2012 at 9:44 PM, Pierre Roudier
<pierre.roudier@gmail.com> wrote:

Thanks Pietro,

Yes that answers the question, but just partly: I was actually
wondering whether it would be possible to get the extracted x and y
coordinates in geographic space (as opposed to the Numpy array space)?

aah sorry, I forgot to say that you can use pixel2coord... see the
example below:

from grass import pygrass
from pygrass.raster import RasterNumpy
from pygrass.region import Region
from pygrass.functions import pixel2coor
elev = RasterNumpy("elevation")
elev.open()
a = elev > 144
reg = Region()
xcoords, ycoords = a.nonzero()
for pixel in zip(list(xcoords), list(ycoords))[:5]:

... print pixel
...
(0, 6)
(0, 7)
(0, 8)
(0, 9)
(0, 10)

for pixel in zip(list(xcoords), list(ycoords))[:5]:

... print pixel2coor(pixel, reg)
...
(228440.0, 630000.0)
(228430.0, 630000.0)
(228420.0, 630000.0)
(228410.0, 630000.0)
(228400.0, 630000.0)

Best regards

Pietro

--
Scientist
Landcare Research, New Zealand

--
Scientist
Landcare Research, New Zealand

Pietro2 · December 18, 2012, 8:46am

Hi Pierre,

On Tue, Dec 18, 2012 at 1:50 AM, Pierre Roudier
<pierre.roudier@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