Opened 15 years ago
Closed 10 years ago
#957 closed defect (fixed)
v.voronoi has extra lines in output
Reported by: | helena | Owned by: | |
---|---|---|---|
Priority: | major | Milestone: | 6.4.4 |
Component: | Vector | Version: | svn-releasebranch64 |
Keywords: | v.voronoi | Cc: | |
CPU: | All | Platform: | All |
Description
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
Attachments (8)
Change History (48)
by , 14 years ago
Attachment: | agac34.zip added |
---|
comment:1 by , 14 years ago
Version: | 6.4.0 RCs → svn-develbranch6 |
---|
by , 14 years ago
Attachment: | Screenshot.png added |
---|
follow-ups: 3 4 comment:2 by , 14 years ago
ping
I can reproduce this bug with 6.4. Any ideas ?
Moritz
comment:3 by , 14 years ago
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
by , 14 years ago
Attachment: | agaccpgr34.zip added |
---|
comment:4 by , 14 years ago
This new file (agaccpgr34.zip) v.voronoi couldn't produce any problem for me. Yalçın
comment:5 by , 14 years ago
Priority: | normal → major |
---|
follow-up: 7 comment:6 by , 13 years ago
Milestone: | 6.4.1 → 6.4.2 |
---|
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
by , 13 years ago
Attachment: | bug_voronoi1.png added |
---|
by , 13 years ago
Attachment: | bug_voronoi2.png added |
---|
follow-up: 8 comment:7 by , 13 years ago
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
follow-ups: 9 19 comment:8 by , 13 years ago
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=boundarysolves 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
by , 13 years ago
Attachment: | bug_voronoi3.png added |
---|
by , 13 years ago
Attachment: | bug_voronoi3.2.png added |
---|
follow-up: 10 comment:9 by , 13 years ago
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
follow-up: 13 comment:10 by , 13 years ago
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
follow-up: 14 comment:11 by , 13 years ago
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
comment:12 by , 13 years ago
Keywords: | v.voronoi added |
---|
comment:13 by , 13 years ago
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.568245610630which 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
follow-up: 15 comment:14 by , 13 years ago
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 thatprintf(%.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
follow-up: 16 comment:15 by , 13 years ago
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
follow-up: 17 comment:16 by , 13 years ago
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
follow-up: 18 comment:17 by , 13 years ago
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!
comment:18 by , 13 years ago
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!
follow-up: 20 comment:19 by , 13 years ago
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
follow-up: 23 comment:20 by , 13 years ago
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
comment:21 by , 12 years ago
Milestone: | 6.4.2 → 6.4.4 |
---|---|
Version: | svn-develbranch6 → svn-releasebranch64 |
follow-up: 24 comment:23 by , 11 years ago
follow-up: 25 comment:24 by , 11 years ago
comment:25 by , 11 years ago
follow-up: 27 comment:26 by , 11 years ago
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
comment:27 by , 11 years ago
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).
follow-ups: 29 30 comment:28 by , 11 years ago
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
comment:29 by , 11 years ago
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
follow-up: 31 comment:30 by , 11 years ago
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.
follow-up: 32 comment:31 by , 11 years ago
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.
follow-ups: 33 36 comment:32 by , 11 years ago
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
follow-up: 34 comment:33 by , 11 years ago
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
follow-up: 35 comment:34 by , 11 years ago
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
comment:35 by , 11 years ago
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.
follow-up: 37 comment:36 by , 11 years ago
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.
follow-up: 38 comment:37 by , 11 years ago
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
comment:38 by , 11 years ago
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).
comment:40 by , 10 years ago
Resolution: | → fixed |
---|---|
Status: | new → closed |
Original problem is solved and no more reports for months. Closing.
I perform voronoi map of 118 sample areas.