Unit 14 - PyGRASS Vector Access

PyGRASS allows to access 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.

Let’s import geodata/osm/counties.gpkg into jena-region GRASS location as described in Unit 03 - Data Management. In this case coordinate reference system (CRS) of input data (EPSG:3035) differs from jena-region CRS definition (EPSG:32632). GRASS will apply transformation of input data into target CRS provided by the current location (Projection match: No). This operation needs to be confirmed by a user, see Fig. 83.

Fig. 83 Confirm reprojection of input data.

Data incosistency can cause significant topological errors when importing vector data into GRASS (as in our case Fig. 84).

Fig. 84 Topological errors detected by import process.

Usually helps to set up a “reasonable” snapping threshold when importing data into GRASS, see Fig. 85.

Fig. 85 Set threshold for snapping before data import in Import settings tab.

Vector data can be treated by Vector class 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 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 ('{0}: {1:.1f}ha'.format(feat.attrs['name'], feat.area() / 10e4))

counties.close()

Example of output:

Hamburg: 31.5ha
Hamburg: 26.6ha
Niedersachsen: 306.4ha
Niedersachsen: 279.2ha
Niedersachsen: 205.7ha
...

Úkol

Modify the script to print overall area for each county.

Writing vector data

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

#!/usr/bin/env python3

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. 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. Topological primitives (areas in this case) are sequentially read by for loop on line 6. Areas are interated by viter() function.
3. For each 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).

#!/usr/bin/env python3

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:

Hamburg                  : 0
Hamburg                  : 1
Niedersachsen            : 1
Niedersachsen            : 1
Niedersachsen            : 1

Úkol

Modify the script to print summary of neighbours for each county.