Unit 14 - PyGRASS Vector Access

PyGRASS allows directly accessing native GRASS raster and vector maps in the sense of Python objects. This unit shows how to deal with GRASS topological vector data by PyGRASS API, see Unit 13 - PyGRASS Raster Access for raster data.

Let’s import geodata/osm/counties.gpkg into jena-region GRASS location.

Vector data

Vector map can be treated by Vector when dealing with basic vector primitives (points, lines, boundaries, centroids) or by VectorTopo which allows accessing topological primitives like areas or isles. See Vector topology section for details.

Example below prints county names and their area size in hectares.

from grass.pygrass.vector import VectorTopo
counties = VectorTopo('counties')
counties.open('r')

for feat in counties.viter('areas'):
    if not feat.attrs:
        continue
    print (u'{0}: {1:.1f}ha'.format(feat.attrs['name'], feat.area() / 10e4))

counties.close()

Writing vector data

In the example below a point vector data with an attribute table will be created.

from grass.pygrass.vector import VectorTopo
from grass.pygrass.vector.geometry import Point

# create the columns definition
cols = [(u'cat',   'INTEGER PRIMARY KEY'),
        (u'name',  'VARCHAR')]
# start new vector with columns definition
new = VectorTopo('pois')
new.open('w', tab_cols=cols, overwrite=True)
# add points
point = Point(681671.15,5644545.63)
new.write(point, ('Jena',))
# commit attributes, otherwise they will be not saved
new.table.conn.commit()
# close the vector
new.close()

Sample script to download: pygrass-write-vector.py

Topology access example

In the following example is presented how to access vector topological primitives directly using PyGRASS. It requires full understanding of GRASS topological model, see Vector topology section in Unit 03 - Data Management.

Sample script below prints for each county number of its neighbours.

  1. Vector map Counties is open on line 3 by VectorTopo and its method open().
  2. Features (areas in this case) are sequentially read by for loop on line 6. Areas are interated by viter() function.
  3. For each feature (ie. county) its boudaries are looped, see line 8. Each boundary has two neighbours (line 9): on the left and right side (-1 for no area).
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
from grass.pygrass.vector import VectorTopo

counties = VectorTopo('counties')
counties.open('r')

for o in counties.viter('areas'):
    neighbours = set()
    for b in o.boundaries():
        for n in b.read_area_ids():
            if n != -1 and n != o.id:
                neighbours.add(n)
    if o.attrs:
        print ('{:25}: {}'.format(o.attrs['name'], len(neighbours)))

counties.close()

Sample script to download: neighbours.py

Possible output:

Baden-Württemberg        : 4
Bayern                   : 6
Saarland                 : 2