How to access vector coordinates in GRASS GIS from python?

Question

Is there a way to iterate over lines and their vertices using python? I'm new to GRASS GIS, and it looks like there are ways to read raster (for numpy), but I don't see a direct way to access coordinates and topology for vectors.

I know that I can access binary data but it does not seem right. I can convert vector feature to db for later use with vector_db_select but it is an extra operation for conversion and I have no data what was what (how can I walk a graph?).

Also when I try to run v.to.db like

v.to.db map=merged_single@PERMANENT layer=1 option=coor units=meters columns=cat

I'm getting

This option requires at least two columns
Finished with error

For instance in ArcGIS I can easily iterate over elements using SearchCursor. I wonder if there is anything similar in GRASS GIS?

Update 2012年06月24日

After reading GRASS GIS manual, it looks like one has to use C library and functions prefixed with Vect_.

Another possible choice probably is to use OGR as it somewhat supports GRASS GIS, but its python bindings are not perfect. It all seems a bit awkward as QGIS so nicely works with all formats from GUI so the inability to access all that from python (which is the official scripting language) seems like a huge shortcoming unless I missed something. Perhaps it will be better in the future.

Answer 1

It is possible to iterate over points/vertices by directly accessing the vector geometry through pygrass.

# Get and plot the coordinates of all of the points in a vector file.
# Then generalize this to lines (this will probably work for boundaries 
# too).
# I'm answering it for both here because (1) this is the way I worked 
# through the problem, and (2) I think it makes intuitive sense because 
# points are simpler and lines are just points connected in an order.
# Written 02 Feb 2015 by Andrew Wickert
# Released to the Public Domain
from grass.pygrass import vector
import numpy as np
from matplotlib import pyplot as plt
############
## POINTS ##
############
data = vector.VectorTopo(NameOfVectorString) # Create a VectorTopo object
data.open('r') # Open this object for reading
# Get the points and append them to the list -- these are point objects.
# These display the coordinates and retain a lot of GRASS GIS functionality.
pointsList = []
for i in range(len(data)):
 pointsList.append(data.read(i+1)) # GRASS vectors seem to be indexed starting at 1
# But for my purposes, it's often nice to just have a numpy array of the actual 
# coordinate values
coords = []
for i in range(len(data)):
 coords.append(data.read(i+1).coords()) # gives a tuple
coords = np.array(coords)
# To see if the points look correct
plt.figure()
plt.plot(coords[:,0], coords[:,1], 'ko')
plt.show()
# And of course you could do better plotting with a tool like Basemap.
# (http://matplotlib.org/basemap/)
###########
## LINES ##
###########
lines = vector.VectorTopo(NameOfVectorString)
lines.open('r')
# On looking at this, each entry in "lines" is a line with multiple points.
# As output, I would like a list of numpy arrays of these points.
def get_vertices(_line):
 """
 Extract vertices in line, in order, and places them in a numpy array
 """
 outline = []
 for vertex in _line:
 outline.append(vertex.coords())
 outline = np.array(outline)
 return outline
LineCoordArrays = [] 
for i in range(len(lines)):
 LineCoordArrays.append( get_vertices(lines[i+1]) )
# Plot them all in whatever color they appear
plt.figure( figsize = (8,8) )
for line in LineCoordArrays:
 plt.plot(line[:,0], line[:,1])
plt.show()

I should mention that this problem has bothered me for quite some time! I have for>1 year used the answer provided by @gene (EDIT: the v.out.ascii solution), and it worked well, but now that I have found this method, I think it is more the "right way". I really appreciate its organization of lines, and while I have not yet tried it with enclosed lines or areas, I expect that it will organize these nicely as well.

Answer 2

For walking a graph you have all the v.net family. To extract the vertices coordinates of lines or polygons (and other columns) one solution is, as always, through v.out.ascii that gives you the xy(z) coordinates(see How one can find the starting and end point of a line in a vector file or how one can connect two line end to end using C code)

test = grass.read_command("v.out.ascii", input="yourvector", format="standard")

The result is given in v.out.ascii. In my case, for example, the variable test is:

ORGANIZATION: 
DIGIT DATE: 
MAP NAME: 
MAP DATE: Sun Jun 24 10:16:35 2012
MAP SCALE: 1
OTHER INFO: 
ZONE: 0
MAP THRESH: 0.000000
VERTI:
L 7 1
206643.21517601 125181.18058576
201007.33432923 121517.8555206
208615.77587567 118699.91687157
199034.77765775 115036.59058769
200725.54321492 111936.8560102
192835.30987687 107428.14794157
192835.30987687 107428.14794157
1 1

It is a (poly)line with one element and 7 vertices (L 7 1) and test is a Python string so

result = test.split("\n") 

and you have a list with the lines of test. If you use regular expressions:

for line in result:
if re.findall(r'^.[0-9]+\.',line):
 print line
 206643.21517601 125181.18058576
 201007.33432923 121517.8555206
 208615.77587567 118699.91687157
 199034.77765775 115036.59058769
 200725.54321492 111936.8560102
 192835.30987687 107428.14794157
 192835.30987687 107428.14794157

and

coords = [] 
for line in result:
 if re.findall(r'^.[0-9]+\.',line): 
 coords.append(line.strip().split(" "))

gives you the solution

[['206643.21517601', '125181.18058576'], ['201007.33432923', '121517.8555206'], ['208615.77587567', '118699.91687157'], ['199034.77765775', '115036.59058769'], ['200725.54321492', '111936.8560102'], ['192835.30987687', '107428.14794157'], ['192835.30987687', '107428.14794157'], ['206643.21517601', '125181.18058576'], ['201007.33432923', '121517.8555206'], ['208615.77587567', '118699.91687157'], ['199034.77765775', '115036.59058769'], ['200725.54321492', '111936.8560102'], ['192835.30987687', '107428.14794157'], ['192835.30987687', '107428.14794157']]

After that, you can create a new point vector file with the xy values and How one can find the starting and end point of a line in a vector file or how one can connect two line end to end using C code

• mlt: I'm accepting this post as an answer as I doubt there are easier approaches. Do you happen to know why v.to.db might fail? While it is not a problem to parse v.out.ascii output, it is definitely neater to work with dbf or something. And I personally going to give OGR python bindings a try first as a more generic approach.
• gene: if you only want to iterate over elements of a layer with v.to.db, see my solution in (osgeo-org.1560.n6.nabble.com/…)
• mlt: Did you paste a wrong link? I don't see that you work with coordinates or the use of v.to.db there.
• gene: It is an example of iteration over a layer with v.db.select. I am working in the Python shell without opening the GRASS application and and I can combine grass.scripts with ogr binding, PyQGIS and modules like Shapely, without problem
• gene: and v.to.db works for points with columns: v.to.db map=pointsmap option=coor col=x,y

Answer 3

If you want to use OGR python bindings, it is easy

from osgeo import ogr
# open grass vector
ds = ogr.Open('/grassdata/geol/test/vector/testline/head')
layer = ds.GetLayer(0)
layer.GetName()
'testline'
 feat = layer.GetFeature(0)
 geom = feature.GetGeometryRef()
 geom.ExportToWkt()
 'LINESTRING (206643.21517600651714 125181.180585758760571,201007.334329231875017 121517.855520597659051,208615.775875667924993 118699.916871574707329,199034.77765774706495 115036.59058768954128,200725.543214922072366 111936.85601020231843,192835.30987687263405 107428.147941565141082,192835.30987687263405 107428.147941565141082)'
 for i in range(geometry.GetPointCount()):
 xy = geometry.GetPoint(i)
 print xy
(206643.21517600652, 125181.18058575876, 0.0)
(201007.33432923188, 121517.85552059766, 0.0)
(208615.77587566792, 118699.91687157471, 0.0)
(199034.77765774706, 115036.59058768954, 0.0)
(200725.54321492207, 111936.85601020232, 0.0)
(192835.30987687263, 107428.14794156514, 0.0)
(192835.30987687263, 107428.14794156514, 0.0)

After that you can use other modules like shapely

from shapely.wkt import loads
line =loads(geometry.ExportToWkt())
list(line.coords) 
[(206643.21517600652, 125181.18058575876), (201007.33432923188, 121517.85552059766), (208615.77587566792, 118699.91687157471), (199034.77765774706, 115036.59058768954), (200725.54321492207, 111936.85601020232), (192835.30987687263, 107428.14794156514), (192835.30987687263, 107428.14794156514)]