v.voronoi has extra lines in output

[helena]

at least since 6.3 v.voronoi produces extra lines in the polygon output, this example run with elev_lid792_randpts has just one, but I have seen worse results:

​http://skagit.meas.ncsu.edu/~helena/grasswork/voronoibug_poly.png

When converted to raster using v.to.rast, the result has holes (nulls), indicating that some polygons do not have values assigned
​http://skagit.meas.ncsu.edu/~helena/grasswork/voronoibug_rast.png

These examples were done with grass65, but I get the same from 64

Helena

[deboyy]

I perform voronoi map of 118 sample areas.

Before perform voronoi for every area set region with the related vector map. But one of the area strange result occured. I check my data and couldn't found any fault or difference to other areas. Is there any bug in the code?

[mlennert]

ping

I can reproduce this bug with 6.4. Any ideas ?

Moritz

[deboyy]

Replying to mlennert:

ping

I can reproduce this bug with 6.4. Any ideas ?

Moritz

Dear Moritz, I attach same areas new shp file with little changes (trees dbh>12 cm). May be it can help to solve problem. I suspect points z coordinate and minimum distance between points.

Yalçın

[deboyy]

This new file (agaccpgr34.zip) v.voronoi couldn't produce any problem for me. Yalçın

[mlennert]

More info on this bug. Using the following data (in file xyz.csv):

x;y;z
0;1;0
1;1;0.8415
2;2;0.9093
3;1;0.1411
4;4;-0.7568
5;3;-0.9589
6;4;-0.2794

and the following commands

v.in.ascii in=xyz.csv out=xyz fs=';' skip=1 x=1 y=2 col="x int, y int, z double precision" --o
g.region vect=xyz
v.voronoi xyz out=xyz_voronoi

I get screenshot bug_voronoi1.png. This seems to show some issue with points at the edges as there are two points located in the same polygon in the bottom left (1 & 2) and top right (6 & 7).

Then I try with a larger region:

g.region n=5 s=0 w=-1 e=7
v.voronoi xyz out=xyz_voronoi2

and I get screenshot bug_voronoi2.png which now seems completely wrong. I've tried this in 6.4.1 and fairly recent versions of 6.5 and trunk (a few days old). The same problem appears in all versions.

Moritz

[mmetz]

Please try trunk r48378. This fixes the extra lines.

The holes reported by Helena are a different problem, a floating point precision error somewhere. Snapping the result, e.g.

v.voronoi in=elev_lid792_randpts out=elev_lid792_randpts_voronoi
v.clean in=elev_lid792_randpts_voronoi out=elev_lid792_randpts_voronoi_snap tool=snap thresh=0.0000001 type=boundary

solves some but not all errors.

Markus M

[mlennert]

Replying to mmetz:

Please try trunk r48378. This fixes the extra lines.

The holes reported by Helena are a different problem, a floating point precision error somewhere. Snapping the result, e.g.

v.voronoi in=elev_lid792_randpts out=elev_lid792_randpts_voronoi
v.clean in=elev_lid792_randpts_voronoi out=elev_lid792_randpts_voronoi_snap tool=snap thresh=0.0000001 type=boundary

solves some but not all errors.

Markus M

For me this solves both my bugs and the one reported by Helena (when going through the extra snap step), i.e. no more extra lines and no holes (after snapping).

However, the following:

g.region vect=precip_30ynormals v.voronoi in=precip_30ynormals out=precip_voronoi

still gives me a polygon containing two points (39 and 82 - see bug_voronoi3.png)

Moritz

[mmetz]

Replying to mlennert:

g.region vect=precip_30ynormals v.voronoi in=precip_30ynormals out=precip_voronoi

still gives me a polygon containing two points (39 and 82 - see bug_voronoi3.png)

The output is correct, no duplicate centroids. I guess in the screenshot you have overlaid precip_30ynormals and not the centroids of precip_voronoi.

The reason why point 82 (and point 112) is missing is that even though g.region vect=precip_30ynormals does set the region extents exactly to the vector extents, the extents are stored with insufficient precision in WIND:

region north: 306221.83019368
vector north: 306221.830193683563

region south: 27606.895351356
vector south: 27606.895351000000

region east: 917004.82916485
vector east: 917004.829164845869

region west: 151768.56824561
vector west: 151768.568245610630

which in this case causes two points to be excluded.

Markus M

[mlennert]

Replying to mmetz:

Replying to mlennert:

g.region vect=precip_30ynormals v.voronoi in=precip_30ynormals out=precip_voronoi

still gives me a polygon containing two points (39 and 82 - see bug_voronoi3.png)

The output is correct, no duplicate centroids. I guess in the screenshot you have overlaid precip_30ynormals and not the centroids of precip_voronoi.

Yes, sorry for not being clear.

The reason why point 82 (and point 112) is missing is that even though g.region vect=precip_30ynormals does set the region extents exactly to the vector extents, the extents are stored with insufficient precision in WIND:

> region north: 306221.83019368
> vector north: 306221.830193683563
> 
> region south: 27606.895351356
> vector south: 27606.895351000000
> 
> region east: 917004.82916485
> vector east: 917004.829164845869
> 
> region west: 151768.56824561
> vector west: 151768.568245610630

which in this case causes two points to be excluded.

Yes, that was it.

g.region vect=precip_30ynormals res=1 -a

does the trick for me.

Should the precision in the WIND file be increased ? Or the value of GRASS_EPSILON to take into account the lesser precision in WIND ?

Moritz

[hamish]

Replying to mmetz:

the extents are stored with insufficient precision in WIND:

> region north: 306221.83019368
> vector north: 306221.830193683563
...

I would suggest that since this is in the 5 nanometer range, and WIND is seemingly stored to %.14g (or does that happen when reading the vector's &window?), that we verify that printf(%.15g) is used when writing the WIND file and then if needed use GRASS_EPSILON to test.

I would suggest ignoring such corner cases in the code, keep using fast >= tests, and recommend in the docs that g.region vect= should be used with the -a and res= options in most cases.

I'm not sure where the %.14g is coming from, as lib/gis/wind_format.c writes using %.15g.

Hamish

[mmetz]

Replying to mlennert:

Replying to mmetz:

Replying to mlennert:

g.region vect=precip_30ynormals v.voronoi in=precip_30ynormals out=precip_voronoi

The reason why point 82 (and point 112) is missing is that even though g.region vect=precip_30ynormals does set the region extents exactly to the vector extents, the extents are stored with insufficient precision in WIND:

> region north: 306221.83019368
> vector north: 306221.830193683563
> 
> region south: 27606.895351356
> vector south: 27606.895351000000
> 
> region east: 917004.82916485
> vector east: 917004.829164845869
> 
> region west: 151768.56824561
> vector west: 151768.568245610630

which in this case causes two points to be excluded.

Yes, that was it.

g.region vect=precip_30ynormals res=1 -a

does the trick for me.

Should the precision in the WIND file be increased ? Or the value of GRASS_EPSILON to take into account the lesser precision in WIND ?

I tried adjusting the bounds with GRASS_EPSILON in g.region, to no avail.

IMHO, this is a more fundamental question. For real projections with units = meter or feet, the current precision is more than enough for raster operations. All vector operations use double precision for coordinates, and here problems start.

1) The reduced precision in the WIND file can cause vector features to fall outside the current region even after g.region vect=myvect.

2) Even if g.region vect=myvect would store the bounds with full precision, from a raster perspective (rows and cols), the vector points falling exactly on the eastern and southern border are outside, not inside the current region. This is because the row index for a region with nrows ranges from 0 to nrows - 1, and (region.north - region.south) / nsres = nrows, but nrows > nrows - 1. A vector point exactly on the southern border would thus be in row number nrows, one row outside the current region. Same for columns and eastern border.

Therefore it may be a good idea to do the adjustment that g.region currently does in

​https://trac.osgeo.org/grass/browser/grass/trunk/general/g.region/main.c#L549

and following lines not only for window.north == window.south etc. but always, to make sure that all vector features are really inside the current region? See also ticket #123.

Markus M

[mmetz]

Replying to hamish:

I would suggest that since this is in the 5 nanometer range, and WIND is seemingly stored to %.14g (or does that happen when reading the vector's &window?), that we verify that printf(%.15g) is used when writing the WIND file and then if needed use GRASS_EPSILON to test.

I would suggest ignoring such corner cases in the code, keep using fast >= tests, and recommend in the docs that g.region vect= should be used with the -a and res= options in most cases.

I'm not sure where the %.14g is coming from, as lib/gis/wind_format.c writes using %.15g.

Not for real projections, at least not for nc_spm_08 with projection: 99 (Lambert Conformal Conic) because full %.15g precision is only used when projection == -1, otherwise %.8f is used. See

​https://trac.osgeo.org/grass/browser/grass/trunk/lib/gis/wind_format.c#L33

Markus M

[hamish]

Replying to hamish:

I'm not sure where the %.14g is coming from, as lib/gis/wind_format.c writes using %.15g.

Replying to mmetz:

Not for real projections, at least not for nc_spm_08 with projection: 99 (Lambert Conformal Conic) because full %.15g precision is only used when projection == -1, otherwise %.8f is used. See ​https://trac.osgeo.org/grass/browser/grass/trunk/lib/gis/wind_format.c#L33

ah, right.. so it's %.8f not %.14g. FWIW I added that -1 for "full" precision* in r37725 as some modules were lying to G_format_*() about not being lat/lon (when they really were) to get decimal output instead of D:M:S, with rather unfortunate results.

[*] (%.15g doesn't quite show full IEEE double, but the max without partial-noise bits)

assuming that projected coords are either in feet or meters, I'm happy enough with 10 nanometer precision. GPS and projection math will take a very long time to catch up with that. Even %.15g won't match what's in the vector's bbox, and again I wonder out loud if this is worth worrying about, or if it should just be declared a feature in the v.voronoi help page.

see also #123.

Hamish

[mmetz]

Replying to hamish:

Even %.15g won't match what's in the vector's bbox, and again I wonder out loud if this is worth worrying about, or if it should just be declared a feature in the v.voronoi help page.

see also #123.

See comment:13 for a suggestion to avoid entries in a number of manuals like "g.region vect=myvect is not enough, you must adjust the region with g.region -a or g.region n=n+x s=s-x e=e+x w=w-x"

Markus M

[hamish]

Replying to mmetz:

See comment:13 for a suggestion

(thanks, I'd missed that)

Therefore it may be a good idea to do the adjustment that g.region currently does in trunk/general/g.region/main.c#L549 and following lines not only for window.north == window.south etc. but always, to make sure that all vector features are really inside the current region?

By "always" you mean just when used with vect=, right? :)

The idea to round out by 1/2 a grid cell in each direction seems ok at first, but I speculate about what the problems might be:

• as res= is typically not set with vect=, what happens if original res is quite large but vector bbox is very small? you get a big region as a result when you wanted one orders of magnitude smaller.

• what if you have a series of vector points broken up by region (lidar 1km x 1km chunk files) and want to patch their bboxes together? would you be better/worse/no worse than now?

perhaps it is better to expand by some other proportional amount? like the lesser of 1/2 a grid cell, or +/- (n-s)*1e-n ? (where n=6..15), or +/- GRASS_EPSILON and change the WIND file (and thus all other G_format_*()) to record to useless/impracticable/wasteful precision? or hardcode the +/- 0.00000001 with a note that it should match the precision in format_double()?

open to suggestions, Hamish

ps- it is wonderful that this longstanding v.voronoi bug is near to being fixed!

[mmetz]

Replying to hamish:

Replying to mmetz:

See comment:13 for a suggestion

(thanks, I'd missed that)

Therefore it may be a good idea to do the adjustment that g.region currently does in trunk/general/g.region/main.c#L549 and following lines not only for window.north == window.south etc. but always, to make sure that all vector features are really inside the current region?

By "always" you mean just when used with vect=, right? :)

Sure :- no absolute always implied here.

The idea to round out by 1/2 a grid cell in each direction seems ok at first, but I speculate about what the problems might be:

• as res= is typically not set with vect=, what happens if original res is quite large but vector bbox is very small? you get a big region as a result when you wanted one orders of magnitude smaller.

Right, but that affects also the current behaviour of g.region

• what if you have a series of vector points broken up by region (lidar 1km x 1km chunk files) and want to patch their bboxes together? would you be better/worse/no worse than now?

The combination of g.region vect=myvect; v.in.region out=myvect_region can exclude features in the current behaviour (%.8f precision). In this case you could use the output of v.info -e plus v.in.ascii to get bboxes and then patch them together (v.info could also get a new output option box which writes the box to a vector).

perhaps it is better to expand by some other proportional amount? like the lesser of 1/2 a grid cell, or +/- (n-s)*1e-n ? (where n=6..15), or +/- GRASS_EPSILON and change the WIND file (and thus all other G_format_*()) to record to useless/impracticable/wasteful precision?

Does it harm to record .15g precision in the WIND file? Internally, it's always stored as double AFAICT.The only harm I can see is implying/suggesting a precision that is not justified by geospatial data. (Unless someone does spatial analysis on electron microscope scans with meters as unit, i.e. non-geo but spatial data).

or hardcode the +/- 0.00000001 with a note that it should match the precision in format_double()?

That one, hardcoding a fixed value matching the minimum precision in the WIND file, makes most sense to me. GRASS_EPSILON has no effect with %.8f since it's 1e-15.

open to suggestions, Hamish

ps- it is wonderful that this longstanding v.voronoi bug is near to being fixed!

[mmetz]

Replying to mlennert:

Replying to mmetz:

Please try trunk r48378. This fixes the extra lines.

The holes reported by Helena are a different problem, a floating point precision error somewhere. Snapping the result, e.g.

> > v.voronoi in=elev_lid792_randpts out=elev_lid792_randpts_voronoi
> > v.clean in=elev_lid792_randpts_voronoi out=elev_lid792_randpts_voronoi_snap tool=snap thresh=0.0000001 type=boundary

solves some but not all errors.

For me this solves both my bugs and the one reported by Helena (when going through the extra snap step), i.e. no more extra lines and no holes (after snapping).

Proper cleaning (a bit more than snapping) is now built into v.voronoi, trunk r48382. The result from v.voronoi in=elev_lid792_randpts out=elev_lid792_randpts_voronoi is now 100% ok.

Markus M

[mmetz]

Replying to mmetz:

Proper cleaning (a bit more than snapping) is now built into v.voronoi, trunk r48382. The result from v.voronoi in=elev_lid792_randpts out=elev_lid792_randpts_voronoi is now 100% ok.

Backported to 6.5 (r48473) and 6.4 (r48474).

The attached test vector agac34.zip is still not working. The bug could be in right_of() in sw_geometry.c.

Markus M

[neteler]

See also #1682

[mmetz]

Replying to mmetz:

The attached test vector agac34.zip is still not working. The bug could be in right_of() in sw_geometry.c.

Fixed in r57268.

[mlennert]

Replying to mmetz:

Replying to mmetz:

The attached test vector agac34.zip is still not working. The bug could be in right_of() in sw_geometry.c.

Fixed in r57268.

AFAICT, this was never backported to grass6. Would it be feasible to do so ?

Moritz

[mmetz]

Replying to mlennert:

Replying to mmetz:

Fixed in r57268.

AFAICT, this was never backported to grass6. Would it be feasible to do so ?

Yes, done in r57819-20.

[hamish]

Hi,

are these empty functions still useful for debug or future work?

​https://trac.osgeo.org/grass/browser/grass/trunk/vector/v.voronoi/sw_output.c#L7

thanks, Hamish

[mmetz]

Replying to hamish:

Hi,

are these empty functions still useful for debug or future work?

​https://trac.osgeo.org/grass/browser/grass/trunk/vector/v.voronoi/sw_output.c#L7

These functions can probably be removed. I guess they come from the original code of Steve J. Fortune (1987).

[mlennert]

In 64release, I'm still getting a subset of the holes reported by Helena in a raster resulting from v.voronoi polygons. Both

v.voronoi input=elev_lid792_randpts@PERMANENT output=voro
v.to.rast input="voro@user1" layer=1 type="point,line,area" output="voro" use="attr" column="value" value=1 rows=4096

and

v.voronoi input=elev_lid792_randpts@PERMANENT output=voro
v.clean in=voro out=voro_snap tool=snap thresh=0.0000001 type=boundary
v.to.rast input="voro_snap@user1" layer=1 type="point,line,area" output="voro_snap" use="attr" column="value" value=1 rows=4096

result in the same holes, although the vector polygons all look ok.

Moritz

[mlennert]

Replying to mlennert:

In 64release, I'm still getting a subset of the holes reported by Helena in a raster resulting from v.voronoi polygons.

Actually I also still see that issue in grass7. The region used in both cases was:

 g.region -p
projection: 99 (Lambert Conformal Conic)
zone:       0
datum:      nad83
ellipsoid:  a=6378137 es=0.006694380022900787
north:      220750
south:      220000
west:       638300
east:       639000
nsres:      1
ewres:      1
rows:       750
cols:       700
cells:      525000

[mmetz]

Replying to mlennert:

In 64release, I'm still getting a subset of the holes reported by Helena in a raster resulting from v.voronoi polygons [...] although the vector polygons all look ok.

That is a bug in v.to.rast, more specifically in G_plot_polygon() in lib/gis/plot.c.

[mmetz]

Replying to mmetz:

Replying to mlennert:

In 64release, I'm still getting a subset of the holes reported by Helena in a raster resulting from v.voronoi polygons [...] although the vector polygons all look ok.

That is a bug in v.to.rast, more specifically in G_plot_polygon() in lib/gis/plot.c.

Please try trunk r57839.

[mlennert]

Replying to mmetz:

Replying to mmetz:

Replying to mlennert:

In 64release, I'm still getting a subset of the holes reported by Helena in a raster resulting from v.voronoi polygons [...] although the vector polygons all look ok.

That is a bug in v.to.rast, more specifically in G_plot_polygon() in lib/gis/plot.c.

Please try trunk r57839.

Solves it for me, thanks !

As this wasn't even a v.voronoi bug, I think this bug report can probably be closed, unless you want to keep it open as a reminder to erase the empty functions. Is there an agreement on erasing them ? If yes, I'm willing to do it.

Moritz

[mlennert]

Replying to mlennert:

Replying to mmetz:

Replying to mmetz:

Replying to mlennert:

In 64release, I'm still getting a subset of the holes reported by Helena in a raster resulting from v.voronoi polygons [...] although the vector polygons all look ok.

That is a bug in v.to.rast, more specifically in G_plot_polygon() in lib/gis/plot.c.

Please try trunk r57839.

Solves it for me, thanks !

Backport to release6 ?

Moritz

[hamish]

Replying to mlennert:

Replying to mlennert:

Replying to mmetz:

Replying to mmetz:

That is a bug in v.to.rast, more specifically in G_plot_polygon() in lib/gis/plot.c.

Please try trunk r57839.

Solves it for me, thanks !

Backport to release6 ?

probably better to test changes to libgis in devbr6 for a little while first.

thanks, Hamish

[mmetz]

Replying to hamish:

Replying to mlennert:

Replying to mlennert:

Replying to mmetz:

Replying to mmetz:

That is a bug in v.to.rast, more specifically in G_plot_polygon() in lib/gis/plot.c.

Please try trunk r57839.

Solves it for me, thanks !

Backport to release6 ?

Backported to 6.5 and 6.4 in r57843,4.

probably better to test changes to libgis in devbr6 for a little while first.

The attempted fix for FPE tolerance has been introduced with GRASS 4.2, but introduced a bug such that sometimes some polygons were not rasterized. This, as the comments say, "weird internal error" underwent quite some testing, I take the liberty to adjust the FPE tolerance in 6.4.svn (BTW, when will we release 6.4.4rc1?) after 16 years of testing.

[mmetz]

Replying to mlennert:

Replying to mmetz:

Replying to mmetz:

Replying to mlennert:

In 64release, I'm still getting a subset of the holes reported by Helena in a raster resulting from v.voronoi polygons [...] although the vector polygons all look ok.

That is a bug in v.to.rast, more specifically in G_plot_polygon() in lib/gis/plot.c.

Please try trunk r57839.

Solves it for me, thanks !

As this wasn't even a v.voronoi bug, I think this bug report can probably be closed, unless you want to keep it open as a reminder to erase the empty functions. Is there an agreement on erasing them ? If yes, I'm willing to do it.

I think these functions can be removed.

[mlennert]

Replying to mmetz:

Replying to mlennert:

Replying to mmetz:

Replying to mmetz:

Replying to mlennert:

In 64release, I'm still getting a subset of the holes reported by Helena in a raster resulting from v.voronoi polygons [...] although the vector polygons all look ok.

That is a bug in v.to.rast, more specifically in G_plot_polygon() in lib/gis/plot.c.

Please try trunk r57839.

Solves it for me, thanks !

As this wasn't even a v.voronoi bug, I think this bug report can probably be closed, unless you want to keep it open as a reminder to erase the empty functions. Is there an agreement on erasing them ? If yes, I'm willing to do it.

I think these functions can be removed.

Done in rXXXXXX for those directly at the spot Hamish mentioned. But I discovered another: ​out_vertex

and two that become empty (except for a conditional test) when the first round is erased: ​out_bisector and ​out_site

I'll attach a diff with all the changes erasing these functions entails. Everything seems to work here, but before submitting this, I would like someone else to have a look and test.

Moritz

This is actually called from

[mlennert]

diff for removing empty functions

[mmetz]

Replying to mlennert: [removing unused functions:]

I'll attach a diff with all the changes erasing these functions entails. Everything seems to work here, but before submitting this, I would like someone else to have a look and test.

I have completely removed all unused functions in r57953 (they were used in the original non-GRASS version producing text output).

[neteler]

Am I right that the ticket can be closed?

[mlennert]

Original problem is solved and no more reports for months. Closing.