[GRASS-dev] [GRASS GIS] #3001: r.in.xyz -s scan mode: one row/col lost?

GDAL · April 25, 2016, 11:05pm

#3001: r.in.xyz -s scan mode: one row/col lost?
----------------------+-------------------------
Reporter: neteler | Owner: grass-dev@…
Type: defect | Status: new
Priority: normal | Milestone: 7.0.4
Component: Raster | Version: svn-trunk
Keywords: r.in.xyz | CPU: All
Platform: Linux |
----------------------+-------------------------
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
}}}

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/3001>
GRASS GIS <https://grass.osgeo.org>

GDAL · April 26, 2016, 1:41am

#3001: r.in.xyz -s scan mode: one row/col lost?
----------------------+-------------------------
  Reporter: neteler | Owner: grass-dev@…
      Type: defect | Status: new
  Priority: normal | Milestone: 7.0.4
Component: Raster | Version: svn-trunk
Resolution: | Keywords: r.in.xyz
       CPU: All | Platform: Linux
----------------------+-------------------------

Comment (by 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:

{{{
#!c
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.

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/3001#comment:1>
GRASS GIS <https://grass.osgeo.org>

GDAL · April 26, 2016, 1:44am

#3001: r.in.xyz -s scan mode: one row/col lost?
----------------------+-------------------------
  Reporter: neteler | Owner: grass-dev@…
      Type: defect | Status: new
  Priority: normal | Milestone: 7.0.4
Component: Raster | Version: svn-trunk
Resolution: | Keywords: r.in.xyz
       CPU: All | Platform: Linux
----------------------+-------------------------
Changes (by wenzeslaus):

* Attachment "r_in_xyz_bounds.png" added.

Result of the original test case (NW corner)

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/3001>
GRASS GIS <https://grass.osgeo.org>

GDAL · April 27, 2016, 8:22am

#3001: r.in.xyz -s scan mode: one row/col lost?
----------------------+-------------------------
  Reporter: neteler | Owner: grass-dev@…
      Type: defect | Status: new
  Priority: normal | Milestone: 7.0.4
Component: Raster | Version: svn-trunk
Resolution: | Keywords: r.in.xyz
       CPU: All | Platform: Linux
----------------------+-------------------------

Comment (by neteler):

Replying to [comment:1 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?

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/3001#comment:2>
GRASS GIS <https://grass.osgeo.org>

GDAL · April 27, 2016, 1:12pm

#3001: r.in.xyz -s scan mode: one row/col lost?
----------------------+-------------------------
  Reporter: neteler | Owner: grass-dev@…
      Type: defect | Status: new
  Priority: normal | Milestone: 7.0.4
Component: Raster | Version: svn-trunk
Resolution: | Keywords: r.in.xyz
       CPU: All | Platform: Linux
----------------------+-------------------------

Comment (by wenzeslaus):

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

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/3001#comment:3>
GRASS GIS <https://grass.osgeo.org>

GDAL · April 29, 2016, 8:54am

#3001: r.in.xyz -s scan mode: one row/col lost?
----------------------+-------------------------
  Reporter: neteler | Owner: grass-dev@…
      Type: defect | Status: new
  Priority: normal | Milestone: 7.0.4
Component: Raster | Version: svn-trunk
Resolution: | Keywords: r.in.xyz
       CPU: All | Platform: Linux
----------------------+-------------------------

Comment (by neteler):

Replying to [comment:3 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

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/3001#comment:4>
GRASS GIS <https://grass.osgeo.org>

GDAL · April 27, 2019, 11:03pm

#3001: r.in.xyz -s scan mode: one row/col lost?
----------------------+-------------------------
  Reporter: neteler | Owner: grass-dev@…
      Type: defect | Status: new
  Priority: normal | Milestone: 7.6.2
Component: Raster | Version: svn-trunk
Resolution: | Keywords: r.in.xyz
       CPU: All | Platform: Linux
----------------------+-------------------------
Changes (by martinl):

* milestone: 7.0.7 => 7.6.2

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/3001#comment:8>
GRASS GIS <https://grass.osgeo.org>