[GRASS-user] Pygrass vector package: write geometric feature in a specific layer

lrntct · December 2, 2014, 3:41pm

Hello all,

I’m writing a GRASS module with pygrass to import data from a text file into a new vector map.
Those data have topological relationships: they are nodes and links of a sewer network. But those objects may have different types of attributes. For example a link could be a pipe or a pump.

Therefore I think that the use of vector layers is the right way to store those data.

With pygrass I managed to set DB links and import the attributes in different tables.
But the only way I’ve found to change the layer of geometric features is with .open(mode=‘w’, layer=X). This works well for one layer, but I’ve to write different objects in different layers.

I tried to close the map and re-open it with another layer, but if I pass mode=‘rw’, the map refuse to open because the layer not yet exist, and with mode=‘w’, it overwrite the existing map.

I know there is a way to create a new map for each layer and then use existing GRASS modules to “merge” them into a single map, but I feel that this is overly complex.

So my question is: Is there a way to decide which layer to write new objects, and change this layer during the import process?

I’m using grass70-svn with the following version of libgis:
libgis Revision: 62395
libgis Date: 2014-10-26 18:17:27

Thanks for you help.

Regards,

Laurent

Pietro2 · December 3, 2014, 6:28am

Hi Laurent,

On Tue, Dec 2, 2014 at 4:41 PM, Laurent C. <lrntct@gmail.com> wrote:

With pygrass I managed to set DB links and import the attributes in
different tables.
But the only way I've found to change the layer of geometric features is
with .open(mode='w', layer=X). This works well for one layer, but I've to
write different objects in different layers.
I tried to close the map and re-open it with another layer, but if I pass
mode='rw', the map refuse to open because the layer not yet exist, and with
mode='w', it overwrite the existing map.

I didn't have the need to do this, so far, so probably I should
rewrite the write method to make it easier.
The first time that you open the map, you should create links and
table and save the link between the vector map and these tables, so
something like (note: this is pseudo-code, not tested! so use it as a
draft idea):

{{{
with VectorTopo('mymap', mode='w', overwrite=True) as vct:
    # create links
    lpipes = Link(layer=1, name='pipes', table=vct.name + 'pipes', key='cat')
    lpumps = Link(layer=2, name='pumps', table=vct.name + 'pumps', key='cat')

    # add the links to the vector map
    if lpipes not in vct.dblinks:
        vct.dblinks.add(lpipes)
    if lpumps not in vct.dblinks:
        vct.dblinks.add(lpumps)

    # create tables removing them if already created
    tpipes = lpipes.table()
    if tpipes.exist():
        tpipes.drop(force=True)
    tpipes.create([('cat', 'PRIMARY KEY'), ('col1', 'VARCHAR(8)'), etc...])
    tpumps = lpumps.table()
    if tpumps.exist():
        tpumps.drop(force=True)
    tpumps.create([('cat', 'PRIMARY KEY'), ('col1', 'VARCHAR(8)'), etc...])

    # write your pipes
    vct.table = tpipes
    for pipe, pipeattr in zip(pipes, pipesattrs):
        vct.write(pipe, attrs=pipeattr)

    # write the pumps
    vct.table = tpumps
    for pump, pumpattr in zip(pumps, pumpattrs):
        # you have to set the category for the second layer
        cats = Cats(pump.c_cats)
        cats.reset()
        cats.set(vct.nlines + 1, lpumps.layer)
        # finally write the pump
        vct.write(pump, attrs=pumpattr, set_cats=False)

}}}

Markus_Metz · December 3, 2014, 7:25am

On Wed, Dec 3, 2014 at 7:28 AM, Pietro <peter.zamb@gmail.com> wrote:

Hi Laurent,

On Tue, Dec 2, 2014 at 4:41 PM, Laurent C. <lrntct@gmail.com> wrote:

With pygrass I managed to set DB links and import the attributes in
different tables.
But the only way I've found to change the layer of geometric features is
with .open(mode='w', layer=X). This works well for one layer, but I've to
write different objects in different layers.
I tried to close the map and re-open it with another layer, but if I pass
mode='rw', the map refuse to open because the layer not yet exist, and with
mode='w', it overwrite the existing map.

I didn't have the need to do this, so far, so probably I should
rewrite the write method to make it easier.
The first time that you open the map, you should create links and
table and save the link between the vector map and these tables, so
something like (note: this is pseudo-code, not tested! so use it as a
draft idea):

How can you write a single feature with multiple links to different
layers? Your example assumes one link per feature. GRASS supports M:N
mapping of vector features to attributes: several features can have
the same category, thus the same attirbutes, and one feature can have
several categories, all in one layer or distributed across several
layers.

Markus M

{{{
with VectorTopo('mymap', mode='w', overwrite=True) as vct:
    # create links
    lpipes = Link(layer=1, name='pipes', table=vct.name + 'pipes', key='cat')
    lpumps = Link(layer=2, name='pumps', table=vct.name + 'pumps', key='cat')

    # add the links to the vector map
    if lpipes not in vct.dblinks:
        vct.dblinks.add(lpipes)
    if lpumps not in vct.dblinks:
        vct.dblinks.add(lpumps)

    # create tables removing them if already created
    tpipes = lpipes.table()
    if tpipes.exist():
        tpipes.drop(force=True)
    tpipes.create([('cat', 'PRIMARY KEY'), ('col1', 'VARCHAR(8)'), etc...])
    tpumps = lpumps.table()
    if tpumps.exist():
        tpumps.drop(force=True)
    tpumps.create([('cat', 'PRIMARY KEY'), ('col1', 'VARCHAR(8)'), etc...])

    # write your pipes
    vct.table = tpipes
    for pipe, pipeattr in zip(pipes, pipesattrs):
        vct.write(pipe, attrs=pipeattr)

    # write the pumps
    vct.table = tpumps
    for pump, pumpattr in zip(pumps, pumpattrs):
        # you have to set the category for the second layer
        cats = Cats(pump.c_cats)
        cats.reset()
        cats.set(vct.nlines + 1, lpumps.layer)
        # finally write the pump
        vct.write(pump, attrs=pumpattr, set_cats=False)

}}}
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

lrntct · December 3, 2014, 9:14pm

2014-12-03 0:28 GMT-06:00 Pietro <peter.zamb@gmail.com>:

you have to set the category for the second layer

cats = Cats(pump.c_cats)
cats.reset()
cats.set(vct.nlines + 1, lpumps.layer)

finally write the pump

vct.write(pump, attrs=pumpattr, set_cats=False)

Hi Pietro,

Thanks for your answer, and for your work on pygrass.

It’s the use of Cats that I was missing. However when I tried that :

{{{
cats = Cats(new_conduit.c_cats)
cats.reset()
cats.set(new_vect_map.nlines + 1, lconduit.layer)
new_vect_map.write(new_conduit, attrs=dbline_content, set_cats=False)
}}}

I got the following error :

"cats.set(new_vect_map.nlines + 1, lconduit.layer)
AttributeError: ‘VectorTopo’ object has no attribute ‘nlines’ "

What ‘nlines’ is supposed to return?

I worked around by using a int variable incremented at each object. It worked fine.

Thanks again !

Regards,

Laurent

···

2014-12-03 0:28 GMT-06:00 Pietro <peter.zamb@gmail.com>:

Hi Laurent,

On Tue, Dec 2, 2014 at 4:41 PM, Laurent C. <lrntct@gmail.com> wrote:

With pygrass I managed to set DB links and import the attributes in
different tables.
But the only way I’ve found to change the layer of geometric features is
with .open(mode=‘w’, layer=X). This works well for one layer, but I’ve to
write different objects in different layers.
I tried to close the map and re-open it with another layer, but if I pass
mode=‘rw’, the map refuse to open because the layer not yet exist, and with
mode=‘w’, it overwrite the existing map.

I didn’t have the need to do this, so far, so probably I should
rewrite the write method to make it easier.
The first time that you open the map, you should create links and
table and save the link between the vector map and these tables, so
something like (note: this is pseudo-code, not tested! so use it as a
draft idea):

{{{
with VectorTopo(‘mymap’, mode=‘w’, overwrite=True) as vct:

create links

lpipes = Link(layer=1, name=‘pipes’, table=vct.name + ‘pipes’, key=‘cat’)
lpumps = Link(layer=2, name=‘pumps’, table=vct.name + ‘pumps’, key=‘cat’)

add the links to the vector map

if lpipes not in vct.dblinks:
vct.dblinks.add(lpipes)
if lpumps not in vct.dblinks:
vct.dblinks.add(lpumps)

create tables removing them if already created

tpipes = lpipes.table()
if tpipes.exist():
tpipes.drop(force=True)
tpipes.create([(‘cat’, ‘PRIMARY KEY’), (‘col1’, ‘VARCHAR(8)’), etc…])
tpumps = lpumps.table()
if tpumps.exist():
tpumps.drop(force=True)
tpumps.create([(‘cat’, ‘PRIMARY KEY’), (‘col1’, ‘VARCHAR(8)’), etc…])

write your pipes

vct.table = tpipes
for pipe, pipeattr in zip(pipes, pipesattrs):
vct.write(pipe, attrs=pipeattr)

write the pumps

vct.table = tpumps
for pump, pumpattr in zip(pumps, pumpattrs):

you have to set the category for the second layer

cats = Cats(pump.c_cats)
cats.reset()
cats.set(vct.nlines + 1, lpumps.layer)

finally write the pump

vct.write(pump, attrs=pumpattr, set_cats=False)

}}}

Pietro2 · December 10, 2014, 9:12am

Hi Laurent,

sorry for the late reply.

On Wed, Dec 3, 2014 at 10:14 PM, Laurent C. <lrntct@gmail.com> wrote:

2014-12-03 0:28 GMT-06:00 Pietro <peter.zamb@gmail.com>:

        # you have to set the category for the second layer
        cats = Cats(pump.c_cats)
        cats.reset()
        cats.set(vct.nlines + 1, lpumps.layer)
        # finally write the pump
        vct.write(pump, attrs=pumpattr, set_cats=False)

Hi Pietro,

Thanks for your answer, and for your work on pygrass.
It's the use of Cats that I was missing. However when I tried that :
{{{
cats = Cats(new_conduit.c_cats)
cats.reset()
cats.set(new_vect_map.nlines + 1, lconduit.layer)
new_vect_map.write(new_conduit, attrs=dbline_content, set_cats=False)
}}}

I got the following error :
"cats.set(new_vect_map.nlines + 1, lconduit.layer)
AttributeError: 'VectorTopo' object has no attribute 'nlines' "

the attribute name was 'n_lines' sorry for the misspelling! and it is
the counter used when writing the geometry feature, basically it does
exactly the thing that you have done manually with the new counter.
However I should change the behavior of the write method to make it
more flexible.

If you have any cleaner ideas on how you would like to interact and
write geometry features may be we can implement it.

Best regards

Pietro

lrntct · July 28, 2025, 7:03pm

Hi,

I hope you don’t mind me digging up this topic, it’s the most recent one I’ve found
I am encountering a strange issue with grass 8.4.1.
I am writing vector geometries—nodes (Points) and link (Line)—and their respective attributes in two different layers and DB tables using the solution from Pietro:

def write_vector_geometry(self, vector_map, geom, cat_num, map_layer):
    """Write geometry in the adequate layer"""
    cats = Cats(geom.c_cats)
    cats.reset()
    cats.set(cat_num, map_layer)
    # write geometry
    vector_map.write(geom, cat_num)

(the attributes are written in a different step)

v.db.connect tells me the layers are correctly connected to each DB table:

v.db.connect -p map=nc_itzi_tutorial_drainage_0003@PERMANENT                    
Vector map <nc_itzi_tutorial_drainage_0003@PERMANENT> is connected by:
layer <1/node> table <nc_itzi_tutorial_drainage_0003_node> in database </home/laurent/data/grassdata/nc_spm_08_grass7/PERMANENT/sqlite/sqlite.db> through driver <sqlite> with key <cat>
layer <2/link> table <nc_itzi_tutorial_drainage_0003_link> in database </home/laurent/data/grassdata/nc_spm_08_grass7/PERMANENT/sqlite/sqlite.db> through driver <sqlite> with key <cat>

Both DB tables are correctly populated with the attributes, and the correct Cat number.

v.db.select shows that the attributes for the link with cat 3 are in the correct DB table:

v.db.select map=nc_itzi_tutorial_drainage_0003@PERMANENT layer=2 format=plain   
cat|link_id|link_type|flow|depth|volume|inlet_offset|outlet_offset|froude
3|C0|conduit|0.000313850743887013|0.00335553119287515|0.0397628607623663|0|0|0.86693500707703

The geometries (in this case, 2 nodes and 1 link) are showing up in the map. But, in grass GUI, the line (cat 3) shows up in layer 1 instead of layer 2, and no attribute is shown in the GUI.

v.category shows this:

v.category input=nc_itzi_tutorial_drainage_0003@PERMANENT option=report         
Layer/table: 1/nc_itzi_tutorial_drainage_0003_node
type       count        min        max
point          2          1          2
line           1          3          3
boundary       0          0          0
centroid       0          0          0
area           0          0          0
face           0          0          0
kernel         0          0          0
all            3          1          3

Doing this:
v.category input=nc_itzi_tutorial_drainage_0003@PERMANENT output=nc_itzi_tutorial_drainage_0003_test@PERMANENT option=chlayer layer=1,2 ids=3

Seems to fix the issue:

v.category input=nc_itzi_tutorial_drainage_0003_test@PERMANENT option=report    
Layer/table: 1/nc_itzi_tutorial_drainage_0003_test_node
type       count        min        max
point          2          1          2
line           0          0          0
boundary       0          0          0
centroid       0          0          0
area           0          0          0
face           0          0          0
kernel         0          0          0
all            2          1          2
Layer/table: 2/nc_itzi_tutorial_drainage_0003_test_link
type       count        min        max
point          0          0          0
line           1          3          3
boundary       0          0          0
centroid       0          0          0
area           0          0          0
face           0          0          0
kernel         0          0          0
all            1          3          3

So, I understand the underlying issue is that when writing the geometric features, the layer information is not properly applied. How could I fix this and make sure the geometries are attached to the correct layer?

lrntct · July 29, 2025, 2:46am

I found a solution to my issue.
While looking at the GRASS source code, I noticed that when writing a new feature to a vector map, the VectorTopo object was using the value of its own “layer” attribute.
So I decided to set this attribute before writing the relevant features. No need to instantiate a Cats object anymore.

Here is a simplified version of my solution:

with VectorTopo(map_name, mode="w", overwrite=self.overwrite) as vector_map:
    # Write in the correct layer
    map_layer, _ = dblinks["node"]
    vector_map.layer = map_layer  # <== here
    for node in drainage_data.nodes:
        point = Point(*node.coordinates)
        # The write function of the vector map sets the layer to the one we set earlier
        vector_map.write(point, cat=cat_num)
        cat_num += 1

    # Set the vector map to the correct layer
    map_layer, _ = dblinks["link"]
    vector_map.layer = map_layer  # <== here
    for link in drainage_data.links:
        line_object = Line(link.vertices)
        # The write function of the vector map sets the layer to the one we set earlier
        vector_map.write(line_object, cat=cat_num)
        cat_num += 1