r.in.xyz: import bug when using scanned extent

[neteler]

When using r.in.xyz in script style, points may get lost during import:

Import script:

# for i in */*.txt
for i in d325095655/p325099657.txt
do
 REGION=`r.in.xyz $i out=dummy fs=space -i --o -sg | cut -d' ' -f1-4`
 g.region $REGION res=1 -p
 r.in.xyz $i out=`basename $i` fs=space -i --o
done

Result:

GRASS 6.3.0svn (pat):~/data/lidar_PAT_raw/raw > sh import_all.sh
projection: 1 (UTM)
zone:       32
datum:      wgs84
ellipsoid:  wgs84
north:      5099457.35
south:      5099071.19
west:       657729.49
east:       657953.44
nsres:      1.00041451
ewres:      0.99977679
rows:       386
cols:       224
cells:      86464
Scanning data ...
Writing to map ...
 100%
r.in.xyz complete. 3 points found in region.

But:

cat d325095655/p325099657.txt
657953.44 5099071.19 542.34
657753.59 5099457.35 327.69
657846.39 5099356.02 736.02
657729.49 5099357.89 585.16

-> one line get's lost.

If I add the "-a" flag to g.region above, it works.

I suspect that a test against GRASS_EPSILON is needed somewhere.

Markus

[hamish]

I will try and have a look at this over the weekend.

Hamish

[neteler]

This is what gets imported (re-imported with g.region res=0.01 -p):

r.stats -1ng p325099657.txt
657753.585 5099457.345 327.6900024414
657729.495 5099357.885 585.1599731445
657846.395 5099356.015 736.0200195312

one point is missing...

[hamish]

Your import_all.sh script should be using -a with res=:

-  g.region $REGION res=1 -p
+  g.region $REGION -a res=1 -p

g.region -a should expand outwards, so no points will be lost.

I will work on a r.in.xyz.auto addon script to scan, set the region, and import all in one step. The problem with that is that res= is unknowable and requires supervision. (I will say more about that in trac bug #37)

Hamish

[hamish]

On grass-dev, Glynn answered:

-> one line get's lost.

Which one? If it's the first line (the southernmost point), that's probably due to this (raster/r.in.xyz/main.c:588):

if (y <= pass_south

y > pass_north) {

The problem with simply changing the test to < is that, for multiple passes, any points on the boundary would be counted twice. The test should be < for the last pass and <= for other passes.

Even for floating-point, copying a value and testing a value for equality with itself should work as expected. Rounding errors don't appear out of nowhere, and neither of those operations should introduce rounding errors.

If GRASS is introducing rounding errors, the first step should be to try to eliminate them (e.g. ensuring that decimal representations have sufficient digits), rather than working around such errors by adding tolerances.

Yes, see this comment in r.in.xyz/main.c ("GE" refers to GRASS_EPSILON)

/* The range should be [0,cols-1]. We use (int) to round down,
     but if the point exactly on eastern edge arr_col will be /just/
     on the max edge .0000000 and end up on the next row.
     We could make above bounds check "if(x>=region.east) continue;"
     But instead we go to all sorts of trouble so that not one single
     data point is lost. GE is too small to catch them all.
     We don't try to make y happy as percent segmenting will make some
     points happen twice that way; so instead we use the y<= test above.
     */

and in description.html:

<h4>Filtering</h4>

Points falling outside the current region will be skipped. This includes
points falling <em>exactly</em> on the southern region bound.
(to capture those adjust the region with "<tt>g.region s=s-0.000001</tt>";
see <em>g.region</em>)

Hamish

[neteler]

Hamish,

if you look at the original report:

I use the scan option to get the region. Using that I lose points. Looks like a precision bug to me (my user hat on). As you can see in my original report, I also used -a for g.region which doesn't make much sense to me if we provide the scan option for creating the g.region parameter string.

Dev hat on I hope that we find a solution for the situation described in "Filtering".

thanks, Markus

[hamish]

I'll have a look at r.in.xyz to see if we can change the < to <= if we are on the last pass block without harming module speed much.

If we do it for all passes there will be duplicated cells for points which fall exactly on pass block borders.

Hamish

[mmetz]

Replying to hamish:

I'll have a look at r.in.xyz to see if we can change the < to <= if we are on the last pass block without harming module speed much.

If we do it for all passes there will be duplicated cells for points which fall exactly on pass block borders.

From a mathematical perspective (is this user hat or dev hat?), the current behaviour is correct.

The user has to decide on a resolution anyway when setting the import region, and in order to get nice clean values, alignment to the resolution is recommended, i.e. have exactly 1 instead of e.g. 0.98990209. In scan mode, the module could print a recommended resolution guestimated from the extends and number of points (works only for evenly spaced points), which can optionally be rounded to integer values.

Since the region settings including resolution and alignment, optionally to the geometry of an existing raster map, are such a core concept of GRASS raster processing, it may not be asked too much to manually set the region including resolution first before importing, and I am inclined to close the ticket as invalid.

Markus M

[neteler]

Closing the ticket as invalid as suggested.