Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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.

share|improve this question
add comment

2 Answers

up vote 1 down vote accepted

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

share|improve this answer
    
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. –  mlt Jun 24 '12 at 23:53
    
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/…) –  gene Jun 25 '12 at 7:30
    
Did you paste a wrong link? I don't see that you work with coordinates or the use of v.to.db there. –  mlt Jun 25 '12 at 7:35
    
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 Jun 25 '12 at 7:39
    
and v.to.db works for points with columns: v.to.db map=pointsmap option=coor col=x,y –  gene Jun 25 '12 at 7:46
add comment

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)]
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.