[Feat] Output ids of points causing fatal error v.net.steiner

Rtutorials commented 3 years ago

Is your feature request related to a problem? Please describe. When the v.net.steiner algorithm fails and says two terminals cannot be connected, the user does not know which points are causing the failure. Nodes are given, but that is script generated and there is no way of determining where those are

Describe the solution you'd like The script prints the feature ids of the points or line segments causing the error

Describe alternatives you've considered showing coordinates

Additional context One such instance could cause the whole thing to fail, see https://gis.stackexchange.com/questions/267624/qgis-v-net-steiner-not-correctly-generated-error

Rtutorials commented 3 years ago

Here is a tiny reproducible example, run v.net.steiner and there is no way of knowing which of the two points failed. Unzip each folder, drag both sqlite files into QGIS, and run v.net.steiner with default settings on the point (examplepoints) and roads (exampleroads) layers

exampleroad.sqlite.zip examplepoints.sqlite.zip

veroandreo commented 3 years ago

Here is a tiny reproducible example, run v.net.steiner and there is no way of knowing which of the two points failed

exampleroad.sqlite.zip examplepoints.sqlite.zip

In general, a really reproducible example includes the commands that need to be run in order to observe the reported behavior/error. This increases the chances to get feedback. Would you mind also providing such commands, please? If possible to be executed in North Carolina sample location.

Rtutorials commented 3 years ago

Edited to add instructions

veroandreo commented 3 years ago

So you are using v.net.steiner from QGIS. Have you tried directly wthin GRASS?

I have no idea of network commands, but I ran v.net.steiner in QGIS 3.18 (using Processing toolbox) and default parameters as suggested. This is the log in order to illustrate the reported problem:

QGIS version: 3.18.2-Zürich
Qt version: 5.15.2
GDAL version: 3.1.4
GEOS version: 3.8.1-CAPI-1.13.3
PROJ version: Rel. 6.3.2, May 1st, 2020
PDAL version: 2.2.0 (git-version: Release)
Processing algorithm…
Algorithm 'v.net.steiner' starting…
Input parameters:
{ '-g' : False, 'GRASS_MIN_AREA_PARAMETER' : 0.0001, 'GRASS_OUTPUT_TYPE_PARAMETER' : 0, 'GRASS_REGION_PARAMETER' : None, 'GRASS_SNAP_TOLERANCE_PARAMETER' : -1, 'GRASS_VECTOR_DSCO' : '', 'GRASS_VECTOR_EXPORT_NOCAT' : False, 'GRASS_VECTOR_LCO' : '', 'acolumn' : '', 'arc_type' : [0,1], 'input' : '/home/veroandreo/Downloads/exampleroad.sqlite|layername=exampleroad', 'npoints' : -1, 'output' : 'TEMPORARY_OUTPUT', 'points' : '/home/veroandreo/Downloads/examplepoints.sqlite|layername=examplepoints', 'terminal_cats' : '1-100000', 'threshold' : 50 }

g.proj -c wkt="/tmp/processing_epZghN/04c835ec74c8437784dc5f9f214845c2/crs.prj"
v.in.ogr min_area=0.0001 snap=-1.0 input="/home/veroandreo/Downloads/exampleroad.sqlite" layer="exampleroad" output="vector_60a2b07d290bd2" --overwrite -o
v.in.ogr min_area=0.0001 snap=-1.0 input="/home/veroandreo/Downloads/examplepoints.sqlite" layer="examplepoints" output="vector_60a2b07d293b23" --overwrite -o
g.region n=3264589.48762953 s=3264496.23383749 e=342015.110712945 w=341737.191647768
v.net -s input=vector_60a2b07d290bd2 points=vector_60a2b07d293b23 output=net60a2b07d29cd04 operation=connect threshold=50.0
v.db.connect -o map=net60a2b07d29cd04 table=vector_60a2b07d293b23 layer=2
v.net.steiner input=net60a2b07d29cd04 arc_type="line,boundary" terminal_cats="1-100000" npoints=-1 output=outputbd663b50926d4a76b1cc9463a71c4827 --overwrite
v.out.ogr type="line" input="outputbd663b50926d4a76b1cc9463a71c4827" output="/tmp/processing_epZghN/d025d4d1375e4c4eafffaf4aa6aa6a2a/output.gpkg" format="GPKG" layer=1 --overwrite
Starting GRASS GIS...
Cleaning up temporary files...
Executing </tmp/processing_epZghN/grassdata/grass_batch_job.sh> ...
Default region was updated to the new projection, but if you have multiple mapsets `g.region -d` should be run in each to update the region from the default
Projection information updated
Over-riding projection check
Check if OGR layer <exampleroad> contains polygons...
0..50..100
Creating attribute table for layer <exampleroad>...
Importing 2 features (OGR layer <exampleroad>)...
0..50..100
-----------------------------------------------------
Building topology for vector map <vector_60a2b07d290bd2@PERMANENT>...
Registering primitives...
Over-riding projection check
Check if OGR layer <examplepoints> contains polygons...
0..33..66..100
Creating attribute table for layer <examplepoints>...
Column name <address ty> renamed to <address_ty>
Column name <geocode ty> renamed to <geocode_ty>
Column name <num of emp> renamed to <num_of_emp>
Column name <employee g> renamed to <employee_g>
Column name <business c> renamed to <business_c>
Column name <monthly te> renamed to <monthly_te>
Column name <monthly ad> renamed to <monthly_ad>
Column name <monthly _1> renamed to <monthly__1>
Column name <monthly _2> renamed to <monthly__2>
Column name <monthly _3> renamed to <monthly__3>
Column name <monthly _4> renamed to <monthly__4>
Column name <monthly _5> renamed to <monthly__5>
Importing 3 features (OGR layer <examplepoints>)...
0..33..66..100
-----------------------------------------------------
Building topology for vector map <vector_60a2b07d293b23@PERMANENT>...
Registering primitives...
Copying features...
50..100
Building topology for vector map <net60a2b07d29cd04@PERMANENT>...
Registering primitives...
Copying attributes...
Building topology for vector map <net60a2b07d29cd04@PERMANENT>...
Registering primitives...
v.net complete. 3 lines (network arcs) written to output.
The table <vector_60a2b07d293b23> is now part of vector map <net60a2b07d29cd04> and may be deleted or overwritten by GRASS modules
Select privileges were granted on the table
Number of terminals: 3
Number of Steiner points set to 1
Building graph...
Registering arcs...
12..25..37..50..62..75..87..100
Flattening the graph...
Graph was built
ERROR: Terminal at node [3] cannot be connected to terminal at node [6]
ERROR: Vector map <outputbd663b50926d4a76b1cc9463a71c4827> not found
Execution of </tmp/processing_epZghN/grassdata/grass_batch_job.sh> finished.
Cleaning up default sqlite database ...
Cleaning up temporary files...
Starting GRASS GIS...
Cleaning up temporary files...
Executing </tmp/processing_epZghN/grassdata/grass_batch_job.sh> ...
ERROR: Vector map <outputbd663b50926d4a76b1cc9463a71c4827> not found
Execution of </tmp/processing_epZghN/grassdata/grass_batch_job.sh> finished.
Cleaning up default sqlite database ...
Cleaning up temporary files...
Execution completed in 1.69 seconds
Results:
{'output': <QgsProcessingOutputLayerDefinition {'sink':TEMPORARY_OUTPUT, 'createOptions': {'fileEncoding': 'System'}}>}

Loading resulting layers
The following layers were not correctly generated.
• /tmp/processing_epZghN/d025d4d1375e4c4eafffaf4aa6aa6a2a/output.gpkg
You can check the 'Log Messages Panel' in QGIS main window to find more information about the execution of the algorithm.

Rtutorials commented 3 years ago

Correct, this is the same error I get. The problem is there is no way to tell which of the three points is causing the failure. Terminal nodes 3 and 6 and just numbers populated by the tool, not the point FIDs. I have never used standalone GRASS.

On Mon, May 17, 2021, 11:14 AM Veronica Andreo @.***> wrote:

So you are using v.net.steiner from QGIS. Have you tried directly wthin GRASS?

I have no idea of network commands, but I ran v.net.steiner in QGIS 3.18 (using Processing toolbox) and default parameters as suggested. This is the log in order to illustrate the reported problem:

QGIS version: 3.18.2-Zürich

Qt version: 5.15.2

GDAL version: 3.1.4

GEOS version: 3.8.1-CAPI-1.13.3

PROJ version: Rel. 6.3.2, May 1st, 2020

PDAL version: 2.2.0 (git-version: Release)

Processing algorithm…

Algorithm 'v.net.steiner' starting…

Input parameters:

{ '-g' : False, 'GRASS_MIN_AREA_PARAMETER' : 0.0001, 'GRASS_OUTPUT_TYPE_PARAMETER' : 0, 'GRASS_REGION_PARAMETER' : None, 'GRASS_SNAP_TOLERANCE_PARAMETER' : -1, 'GRASS_VECTOR_DSCO' : '', 'GRASS_VECTOR_EXPORT_NOCAT' : False, 'GRASS_VECTOR_LCO' : '', 'acolumn' : '', 'arc_type' : [0,1], 'input' : '/home/veroandreo/Downloads/exampleroad.sqlite|layername=exampleroad', 'npoints' : -1, 'output' : 'TEMPORARY_OUTPUT', 'points' : '/home/veroandreo/Downloads/examplepoints.sqlite|layername=examplepoints', 'terminal_cats' : '1-100000', 'threshold' : 50 }

g.proj -c wkt="/tmp/processing_epZghN/04c835ec74c8437784dc5f9f214845c2/crs.prj"

v.in.ogr min_area=0.0001 snap=-1.0 input="/home/veroandreo/Downloads/exampleroad.sqlite" layer="exampleroad" output="vector_60a2b07d290bd2" --overwrite -o

v.in.ogr min_area=0.0001 snap=-1.0 input="/home/veroandreo/Downloads/examplepoints.sqlite" layer="examplepoints" output="vector_60a2b07d293b23" --overwrite -o

g.region n=3264589.48762953 s=3264496.23383749 e=342015.110712945 w=341737.191647768 v.net -s input=vector_60a2b07d290bd2 points=vector_60a2b07d293b23 output=net60a2b07d29cd04 operation=connect threshold=50.0

v.db.connect -o map=net60a2b07d29cd04 table=vector_60a2b07d293b23 layer=2

v.net.steiner input=net60a2b07d29cd04 arc_type="line,boundary" terminal_cats="1-100000" npoints=-1 output=outputbd663b50926d4a76b1cc9463a71c4827 --overwrite

v.out.ogr type="line" input="outputbd663b50926d4a76b1cc9463a71c4827" output="/tmp/processing_epZghN/d025d4d1375e4c4eafffaf4aa6aa6a2a/output.gpkg" format="GPKG" layer=1 --overwrite

Starting GRASS GIS...

Cleaning up temporary files...

Executing </tmp/processing_epZghN/grassdata/grass_batch_job.sh> ...

Default region was updated to the new projection, but if you have multiple mapsets g.region -d should be run in each to update the region from the default

Projection information updated

Over-riding projection check

Check if OGR layer contains polygons...

0..50..100

Creating attribute table for layer ...

Importing 2 features (OGR layer )...

0..50..100

Building topology for vector map vector_60a2b07d290bd2@PERMANENT...

Registering primitives...

Over-riding projection check

Check if OGR layer contains polygons...

0..33..66..100

Creating attribute table for layer ...

Column name
renamed to

Column name renamed to

Column name renamed to

Column name renamed to

Column name renamed to

Column name renamed to

Column name renamed to

Column name renamed to

Column name renamed to

Column name renamed to

Column name renamed to

Column name renamed to

Importing 3 features (OGR layer )...

0..33..66..100

Building topology for vector map vector_60a2b07d293b23@PERMANENT...

Registering primitives...

Copying features...

50..100

Building topology for vector map net60a2b07d29cd04@PERMANENT...

Registering primitives...

Copying attributes...

Building topology for vector map net60a2b07d29cd04@PERMANENT...

Registering primitives... v.net complete. 3 lines (network arcs) written to output.

The table is now part of vector map and may be deleted or overwritten by GRASS modules

Select privileges were granted on the table

Number of terminals: 3

Number of Steiner points set to 1

Building graph...

Registering arcs...

12..25..37..50..62..75..87..100

Flattening the graph...

Graph was built

ERROR: Terminal at node [3] cannot be connected to terminal at node [6]

ERROR: Vector map not found

Execution of </tmp/processing_epZghN/grassdata/grass_batch_job.sh> finished.

Cleaning up default sqlite database ...

Cleaning up temporary files...

Starting GRASS GIS...

Cleaning up temporary files...

Executing </tmp/processing_epZghN/grassdata/grass_batch_job.sh> ...

ERROR: Vector map not found

Execution of </tmp/processing_epZghN/grassdata/grass_batch_job.sh> finished.

Cleaning up default sqlite database ...

Cleaning up temporary files...

Execution completed in 1.69 seconds

Results:

{'output': <QgsProcessingOutputLayerDefinition {'sink':TEMPORARY_OUTPUT, 'createOptions': {'fileEncoding': 'System'}}>}

Loading resulting layers

The following layers were not correctly generated.

• /tmp/processing_epZghN/d025d4d1375e4c4eafffaf4aa6aa6a2a/output.gpkg

You can check the 'Log Messages Panel' in QGIS main window to find more information about the execution of the algorithm.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/OSGeo/grass/issues/1572#issuecomment-842530943, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEXIH6FC266TY6ITQEGIIATTOFMIZANCNFSM44XBWFQQ .

neteler commented 3 years ago

Edited to add instructions

Perhaps our definition of "instructions" is something else but luckily we are used to analyse unknown data :-)

Prepare data

# check unknown data, is it geospatial at all?
ogrinfo -so -al examplepoints.sqlite
ogrinfo -so -al exampleroad.sqlite

# create GRASS GIS location from dataset
grass78 -c examplepoints.sqlite  ~/grassdata/utm17n_test

# import data
v.import examplepoints.sqlite
v.import exampleroad.sqlite

g.region vector=exampleroad -p
g.gui

Network preparation (connect nodes to lines), with snapping enabled (2m threshold):

(as per documentation, a mandatory step!)

v.net -s input=exampleroad points=examplepoints output=network operation=connect threshold=2

Note that nodes are on layer 2, hence one have to enable "show all layers" to see them:

Check attributes

v.db.select network  layer=-1
cat|ogc_fid2|ogc_fid1|ogc_fid0|linearid|fullname|rttyp|mtfcc|layer|path|networkgrp
1|1|2864|905|1106087315821|NE 115th St|M|S1400|tl_2020_12075_roads|/Users/REDACTED/Downloads/FLrds/tl_2020_12075_roads.shp|0
2|2|3160|980|1103698308515|NE 85th Ave|M|S1400|tl_2020_12075_roads|/Users/REDACTED/Downloads/FLrds/tl_2020_12075_roads.shp|0

Network analysis

Using v.net.steiner (manual):

(hint: to show node numbers on the map enable in the vector display [x] Display topology information (nodes, edges))

v.net.steiner network arc_layer=1 node_layer=2 output=steiner terminal_cats=1-2
Number of terminals: 2
Number of Steiner points set to 0
Building graph...
Registering arcs...
 100%
Flattening the graph...
Graph was built
[3] (not reachable) nodes removed from list of Steiner point candidates
MST costs = 14.355000
Number of added Steiner points: 0 (theoretic max is 0).
Steiner tree costs = 14.355000
Building topology for vector map <steiner@PERMANENT>...
Registering primitives...
A Steiner tree with 1 arc has been built

So far, so nice.

neteler commented 3 years ago

Indeed, when trying a different set of connections it fails:

v.net.steiner network arc_layer=1 node_layer=2 output=steiner terminal_cats=1-6
Number of terminals: 3
Number of Steiner points set to 1
Building graph...
Registering arcs...
 100%
Flattening the graph...
Graph was built
ERROR: Terminal at node [3] cannot be connected to terminal at node [6]

This is the error message which you have reported.

The question is: is it the data or the algorithm?

Rtutorials commented 3 years ago

This is the image of the network layer with View all layers set to -1. When it says it can't connect 3 and 6, does that refer to n3 and n6 on the map?

The error itself is likely the data, running a dissolve in QGIS fixes the issue for most but not all of these instances

neteler commented 3 years ago

More detailed output

BTW, one can get more details with "mild" debugging:

# deactivate debugging
g.gisenv set="DEBUG=1"

# run Steiner network analysis
v.net.steiner network arc_layer=1 node_layer=2 output=steiner terminal_cats=1-6 --o
D1/1: G_set_program_name(): v.net.steiner
D1/1: Imput categories:
D1/1: 1 - 6
...
Number of terminals: 3
Number of Steiner points set to 1
D1/1: List of terminal nodes (3):

D1/1: 3
D1/1: 2
D1/1: 6

D1/1: Vect_net_build_graph(): ltype = 6, afield = 1, nfield = 0
D1/1:     afcol = (null), abcol = (null), ncol = (null)
Building graph...
Registering arcs...
 100%
Flattening the graph...
Graph was built
D1/1: Destination node 5 is unreachable from node 3
D1/1: Destination node 6 is unreachable from node 3
D1/1: Destination node 7 is unreachable from node 3
D1/1: Destination node 5 is unreachable from node 2
D1/1: Destination node 6 is unreachable from node 2
D1/1: Destination node 7 is unreachable from node 2
D1/1: Destination node 1 is unreachable from node 6
D1/1: Destination node 2 is unreachable from node 6
D1/1: Destination node 3 is unreachable from node 6
D1/1: Destination node 4 is unreachable from node 6
ERROR: Terminal at node [3] cannot be connected to terminal at node [6]

# deactivate debugging
g.gisenv set="DEBUG=0"

With level 3 is is even more details:

g.gisenv set="DEBUG=3"

v.net.steiner network arc_layer=1 node_layer=2 output=steiner terminal_cats=1-6 --o
...
D1/3: Destination node 1 is unreachable from node 6
D3/3: init costs 6 - > 1 = -2.000000
D3/3: find_shortest_path(): from = 6, to = 2
D1/3: Destination node 2 is unreachable from node 6
D3/3: init costs 6 - > 2 = -2.000000
D3/3: find_shortest_path(): from = 6, to = 3
D1/3: Destination node 3 is unreachable from node 6
D3/3: init costs 6 - > 3 = -2.000000
D3/3: find_shortest_path(): from = 6, to = 4
D1/3: Destination node 4 is unreachable from node 6
D3/3: init costs 6 - > 4 = -2.000000
D3/3: find_shortest_path(): from = 6, to = 5
D3/3: init costs 6 - > 5 = 45.948000
D3/3: find_shortest_path(): from = 6, to = 7
D3/3: init costs 6 - > 7 = 245.108000
ERROR: Terminal at node [3] cannot be connected to terminal at node [6]

# deactivate debugging
g.gisenv set="DEBUG=0"

A Steiner tree from this sample dataset

To me it appears to be a problem of the dataset.

v.net.steiner network arc_layer=1 node_layer=2 output=steiner terminal_cats=1-2 --o 
Number of terminals: 2
Number of Steiner points set to 0
Building graph...
Registering arcs...
 100%
Flattening the graph...
Graph was built
[3] (not reachable) nodes removed from list of Steiner point candidates
MST costs = 14.355000
Number of added Steiner points: 0 (theoretic max is 0).
Steiner tree costs = 14.355000
WARNING: Vector map <steiner> already exists and will be overwritten
Building topology for vector map <steiner@PERMANENT>...
Registering primitives...
A Steiner tree with 1 arc has been built

Steiner Tree found (shown in Orange):

mlennert commented 3 years ago

This is the image of the network layer with View all layers set to -1. When it says it can't connect 3 and 6, does that refer to n3 and n6 on the map?

I believe so.

OSGeo / grass

[Feat] Output ids of points causing fatal error v.net.steiner #1572

More detailed output

A Steiner tree from this sample dataset