Wednesday 16 May 2018

qgis - How to access vector coordinates in GRASS GIS from python?


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



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


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...