r.in.xyz -s scan mode: one row/col lost?

[neteler]

While preparing an example how to import X,Y,Z ASCII DEMs which look like this:

630007.5 228492.5 141.99614
630022.5 228492.5 141.37904
630037.5 228492.5 142.29822
630052.5 228492.5 143.97987
...

for the page

https://grass.osgeo.org/grass70/manuals/r.in.xyz.html#import-of-x,y,string-data

I came across the issue that one row/col seems to be lost in the scan mode of r.in.xyz. If true this may have some relevance for Lidar processing, too.

Test case:

#### GRASS 7.0.4svn (nc_spm_08_grass7):~ > 

# export as X,Y,Z ASCII
g.region raster=elevation -p
projection: 99 (Lambert Conformal Conic)
zone:       0
datum:      nad83
ellipsoid:  a=6378137 es=0.006694380022900787
north:      228500
south:      215000
west:       630000
east:       645000
nsres:      10
ewres:      10
rows:       1350
cols:       1500
cells:      2025000

# export (cell centers)
r.stats -1g elevation > elevation.xyz

head -n 4 elevation.xyz
630005 228495 141.99614
630015 228495 141.27849
630025 228495 141.37904
630035 228495 142.29822

wc -l elevation.xyz
2025000 elevation.xyz

# ... looks ok.

#### import
# scan extent from xyz input
r.in.xyz input=elevation.xyz separator=space -s output=dummy -g
n=228495 s=215005 e=644995 w=630005 b=55.578793 t=156.32986

# set region
g.region n=228495 s=215005 e=644995 w=630005 -p
projection: 99 (Lambert Conformal Conic)
zone:       0
datum:      nad83
ellipsoid:  a=6378137 es=0.006694380022900787
north:      228495
south:      215005
west:       630005
east:       644995
nsres:      10
ewres:      10
rows:       1349
cols:       1499
cells:      2022151

## -->> one col and one row lost?!

r.in.xyz input=elevation.xyz separator=space method=mean output=myelev
...
r.in.xyz complete. 2022151 points found in region.

## -->> incomplete

[wenzeslaus]

The exported points are in a grid and the points are placed in centers of original raster cells. This explains why the region is one cell smaller. With the same resolution and alignment to the input data (i.e. to the points rather than the original raster). This causes the new raster to begin in the middle of the cell from the original raster and end similarly in the middle of the original cell. It is therefore one cell smaller (2 * 0.5 cell) than the original.

During the actual import, the points on some borders (I was not able to figure out which ones) are probably left out due to floating point comparisons. When you make the region 1 meter bigger on each side, all points are imported and the central row and column have 2 cells in them and central cell has 4.

Here is the relevant r.in.xyz code:

if (y <= region.south || y > region.north) {
    G_free_tokens(tokens);
    continue;
}
if (x < region.west || x >= region.east) {
    G_free_tokens(tokens);
    continue;
}
arr_row = (int)((region.north - y) / region.ns_res) - row0;
if (arr_row < 0 || arr_row >= rows) {
    G_free_tokens(tokens);
    continue;
}
arr_col = (int)((x - region.west) / region.ew_res);
/* no check of arr_col */

r.in.lidar behaves the same as r.in.xyz and uses the same code to exclude the points. Tested with:

r.to.vect in=elevation out=elevation_vector type=point -btz
v.out.lidar in=elevation_vector out=elevation.las

v.in.ascii -r (import in region only) gives all the points, the code is:

if ((window.east < easting) ||
    (window.west > easting))
    skip = TRUE;
if ((window.north < northing) ||
    (window.south > northing))
    skip = TRUE;

Changing the <= and >= to < and > imports 2023500 points (more than 2022151 but not 2025000) but that doesn't look like good idea anyway. I'm even worried about not checking arr_col (I ended up checking row, col, depth in r3.in.lidar since checking coordinates was not reliable).

I think the current behavior might miss few points (max 4?) for randomly distributed points when region would be aligned to points rather than resolution (g.region -a). For points in grid as in this case, the user is responsible for dealing with floating point issues and gridded data specifics with the current implementation.

[wenzeslaus]

Result of the original test case (NW corner)

[neteler]

Replying to wenzeslaus: ...

For points in grid as in this case, the user is responsible for dealing with floating point issues and gridded data specifics with the current implementation.

Well, IMHO 99% of the users will not be able to deal with this issue for points in grids. Perhaps we need a simple addon wrapping r.in.xyz to get it right? Or a flag to deal properly with import of points in grids?

[wenzeslaus]

Flag or special module might be possible, but what would they do? Identify the cell size of the grid?

[neteler]

Replying to wenzeslaus:

Flag or special module might be possible, but what would they do? Identify the cell size of the grid?

That would be ideal but probably hard to implement.

For now I have improved the example https://grass.osgeo.org/grass70/manuals/r.in.xyz.html#import-of-x,y,z-ascii-into-dem