Introduction to Vector classes

Details about the GRASS GIS vector architecture can be found in the GRASS GIS 7 Programmer’s Manual: GRASS Vector Library

PyGrass has two classes for vector maps: Vector and VectorTopo. As the names suggest, the Vector class is for vector maps, while VectorTopo opens vector maps with GRASS GIS topology. VectorTopo is an extension of the Vector class, so supports all the Vector class methods, with additions.

Vector

The Vector class is part of the vector module. It is based on the Info class, which provides methods for accessing basic information about the vector map:

>>> from grass.pygrass.vector import Vector
>>> cens = Vector('census')
>>> cens.is_open()
False
>>> cens.mapset
''
>>> cens.exist()
True
>>> cens.mapset
'PERMANENT'
>>> cens.overwrite
False

VectorTopo

The VectorTopo class allows vector maps to be loaded along with an accompanying topology. The VectorTopo class interface has all the same methods as the basic Vector class, but includes many more methods dependent on the topology rules. Consult each class’s documentation for a list of the methods to compare what each class can do. Using the VectorTopo class is just like the Vector class:

>>> from grass.pygrass.vector import VectorTopo
>>> municip = VectorTopo('boundary_municp_sqlite')
>>> municip.is_open()
False
>>> municip.mapset
''
>>> municip.exist()  # check if exist, and if True set mapset
True
>>> municip.mapset
'user1'

Working with Vector Objects

As the VectorTopo class is so similar to the Vector class, the following examples exclusively demonstrate the VectorTopo class.

To begin using a vector map, it must first be opened:

>>> from grass.pygrass.vector import VectorTopo
>>> municip = VectorTopo('boundary_municp_sqlite')
>>> municip.open(mode='r')

The open() method supports a number of option arguments (see the Info documentation for a complete list). In particular, the mode argument can take a a value of r for reading, w for writing, or rw for reading/writing.

The geometry of a vector map can be read sequentially using the next() method. To return to the beginning, use the rewind() method.

>>> municip.next()
Boundary(v_id=1)
>>> municip.next()
Boundary(v_id=2)
>>> municip.next()
Boundary(v_id=3)
>>> municip.rewind()
>>> municip.next()
Boundary(v_id=1)

If a vector map is opened with the mode w or rw, then the user can write new features to the dataset:

Open a new vector map:

>>> new = VectorTopo('newvect')
>>> new.exist()
False

Define the new columns in the attribute table:

>>> cols = [(u'cat',       'INTEGER PRIMARY KEY'),
...         (u'name',      'TEXT')]

Open the vector map in write mode:

>>> new.open('w', tab_name='newvect', tab_cols=cols)

Import the geometry feature class and add two points:

>>> from grass.pygrass.vector.geometry import Point
>>> point0 = Point(636981.336043, 256517.602235)
>>> point1 = Point(637209.083058, 257970.129540)

Write the two points to the map:

>>> new.write(point0, ('pub', ))
>>> new.write(point1, ('resturnat', ))

Commit the db changes:

>>> new.table.conn.commit()
>>> new.table.execute().fetchall()
[(1, u'pub'), (2, u'resturnat')]

Close the vector map:

>>> new.close()
>>> new.exist()
True

Now we can play with the map:

>>> new.open(mode='r')
>>> new.read(1)
Point(636981.336043, 256517.602235)
>>> new.read(2)
Point(637209.083058, 257970.129540)
>>> new.read(1).attrs['name']
u'pub'
>>> new.read(2).attrs['name']
u'resturnat'
>>> new.close()
>>> new.remove()

Note the close() and remove() methods above. The close() method ensure that the files are properly released by the operating system, and will also build the topology (if called by the VectorTopo class). Take caution with the remove() method; it is used here because this example map was temporary. Calling this method will completely delete the map and its data from the file system.

More Features of the VectorTopo Class

See the class documentation for a full list of methods, but here are some examples using VectorTopo methods:

Get the number of primitives:

>>> municip.num_primitive_of('line')
0
>>> municip.num_primitive_of('centroid')
3579
>>> municip.num_primitive_of('boundary')
5128

Get the number of different feature types in the vector map:

>>> municip.number_of("areas")
3579
>>> municip.number_of("islands")
2629
>>> municip.number_of("holes")
0
>>> municip.number_of("lines")
8707
>>> municip.number_of("nodes")
4178
>>> municip.number_of("pizza")  
Traceback (most recent call last):
    ...
ValueError: vtype not supported, use one of: 'areas', ..., 'volumes'
Traceback (most recent call last):
    ...
ValueError: vtype not supported, use one of: 'areas', ..., 'volumes'

Note that the method with raise a ValueError if a non-supported vtype is specified.

Accessing Attribute Tables

The GRASS philosophy stipulates that vector map features are independent from their attributes, and that a vector map’s attribute table(s) should not be loaded unless explicitly specified, in case they are not needed.

Accessing a vector map’s table(s) requires finding any links to tables, then requesting the table from each of the returned links:

>>> from grass.pygrass.vector import VectorTopo
>>> municip = VectorTopo('census')
>>> municip.open(mode='r')
>>> dblinks = DBlinks(municip.c_mapinfo)
>>> dblinks
DBlinks([Link(1, census, sqlite)])
>>> link = DBlinks[0]
Link(1, census, sqlite)
>>> table = link.table()

Here, DBlinks() is a class (DBlinks) that contains all the links of a vector map. Each link is also a class (Link) that contains a specific link’s parameters. The table() method of the link class return the linked table as a table object (Table).

Geometry Classes

The vector package also includes a number of geometry classes, including Area, Boundary, Centroid, Isle, Line, and Point classes. Please consult the geometry module for a complete list of methods for these classes, as there are many. Some basic examples are given below.

Instantiate a Point object that could be 2 or 3D, default parameters are 0:

>>> pnt = Point()
>>> pnt.x
0.0
>>> pnt.y
0.0
>>> pnt.z
>>> pnt.is2D
True
>>> pnt
Point(0.000000, 0.000000)
>>> pnt.z = 0
>>> pnt.is2D
False
>>> pnt
Point(0.000000, 0.000000, 0.000000)
>>> print(pnt)
POINT(0.000000 0.000000 0.000000)

Create a Boundary and calculate its area:

>>> bound = Boundary(points=[(0, 0), (0, 2), (2, 2), (2, 0),
...                          (0, 0)])
>>> bound.area()
4.0

Construct a Line feature and find its bounding box:

>>> line = Line([(0, 0), (1, 1), (2, 0), (1, -1)])
>>> line
Line([Point(0.000000, 0.000000),
      Point(1.000000, 1.000000),
      Point(2.000000, 0.000000),
      Point(1.000000, -1.000000)])
>>>bbox = line.bbox()
>>> bbox
Bbox(1.0, -1.0, 2.0, 0.0)

Buffer a Line feature and find the buffer centroid:

>>> line = Line([(0, 0), (0, 2)])
>>> area = line.buffer(10)
>>> area.boundary
Line([Point(-10.000000, 0.000000),...Point(-10.000000, 0.000000)])
>>> area.centroid
Point(0.000000, 0.000000)

More Examples

Find all areas larger than 10000m2:

>>> big = [area for area in municip.viter('areas')
...        if area.alive() and area.area >= 10000]

The PyGrass vector methods make complex operations rather easy. Notice the viter() method: this returns an iterator object of the vector features, so the user can choose on which vector features to iterate without loading all the features into memory.

We can then sort the areas by size:

>>> from operator import methodcaller as method
>>> big.sort(key = method('area'), reverse = True)  # sort the list
>>> for area in big[:3]:
...     print area, area.area()
Area(3102) 697521857.848
Area(2682) 320224369.66
Area(2552) 298356117.948

Or sort for the number of isles that are contained inside:

>>> big.sort(key = lambda x: x.isles.__len__(), reverse = True)
>>> for area in big[:3]:
...     print area, area.isles.__len__()
...
Area(2682) 68
Area(2753) 45
Area(872) 42

Or can list only the areas containing isles:

>>> area_with_isles = [area for area in big if area.isles]
>>> area_with_isles                                   
[Area(...), ..., Area(...)]

Of course is still possible work only with a specific area, with:

>>> from grass.pygrass.vector.geometry import Area
>>> area = Area(v_id=1859, c_mapinfo=municip.c_mapinfo)
>>> area.area()
39486.05401495844
>>> area.bbox()  # north, south, east, west
Bbox(175711.718494, 175393.514494, 460344.093986, 460115.281986)
>>> area.isles
Isles([])

Now, find an area with an island inside...

>>> area = Area(v_id=2972, c_mapinfo=municip.c_mapinfo)
>>> area.isles                                       
Isles([Isle(1538), Isle(1542), Isle(1543), ..., Isle(2571)])
>>> isle = area.isles[0]
>>> isle.bbox()
Bbox(199947.296494, 199280.969494, 754920.623987, 754351.812986)