[GRASS-dev] [GRASS GIS] #1385: r3.in.ascii man page vs. r3.in.ascii source code

GRASS GRASS developer
1

#1385: r3.in.ascii man page vs. r3.in.ascii source code
-------------------------+--------------------------------------------------
Reporter: hamish | Owner: grass-dev@…
     Type: defect | Status: new
Priority: normal | Milestone: 6.4.2
Component: Docs | Version: svn-develbranch6
Keywords: r3.in.ascii | Platform: Linux
      Cpu: x86-64 |
-------------------------+--------------------------------------------------
Hi,

the man page for `r3.in.ascii` (6.5svn) states:

  ''Note lower-left (SW) corner of the bottom level comes first! This
array format, where EW is preserved but NS is flipped, is some-times known
as "ij" coordinates. This is opposite to r.in.ascii's format, which
places the SW corner at the beginning of the last row of data.''

but then shows data format example like:
{{{
,y1,
,y2,
,y3,
}}}
(although I think I might have added that :slight_smile:

and importing data with NW corner as the first bytes of the file, then
`r3.to.rast` to split into layers seems to successfully display north-up.
So either the man page or the `r3.in.ascii` code has got it backwards..?

I also notice that `r3.out.vtk` exports celldata with the data south-up
within the ascii.vtk file, but I'm not sure if that is part of the VTK
design, or the export module flipping what it doesn't actually have to. ??

----
tangent:
I don't know Paraview well enough to reliably say which way 'should' be up
there. The first thing I notice in Paraview 3.8.0 after loading the data,
apply, Display|Style->representation,Surface + Color->Color by->''layer
name'' is that the top level is on the bottom and the bottom layer is on
the top. again in the display tab setting Transformation->Scale->z to -1
makes it ok, but it gets very dark, as if it is all in shadow. (test data
is cube containing 0/1 boolean values, similar to what r.to.rast3elev
would create from a 2.5D DEM for land/air. I'm still struggling to
visualize that well in NVIZ or paraview.)
??or perhaps I set top: and bottom: backwards in the `r3.in.ascii` header?
`r3.to.rast` shows level 0001 as top layer, level 25 (of 25) as bottom
layer. (data is all below sea level)

Hamish

ps- don't try searching in trac for `r3.in.ascii`, as it treats the `r3`
as svn rev3 which is almost every file in the grass 5 source code! The web
browser does not enjoy it.

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/1385&gt;
GRASS GIS <http://grass.osgeo.org>

2

#1385: r3.in.ascii man page vs. r3.in.ascii source code
-------------------------+--------------------------------------------------
Reporter: hamish | Owner: grass-dev@…
     Type: defect | Status: new
Priority: normal | Milestone: 6.4.2
Component: Docs | Version: svn-develbranch6
Keywords: r3.in.ascii | Platform: Linux
      Cpu: x86-64 |
-------------------------+--------------------------------------------------

Comment(by huhabla):

Hi,
the internal addressing of voxel data in the g3d lib is equal to the
raster library in x,y coordinates -> the counting starts at upper left
corner. The depths are counted from the bottom to the top. All modules
which i have implemented following this approach.
You need to look deep in the code, because this behavior is not explicitly
documented, but this code from the g3d lib makes it clear:

{{{
void
G3d_getRegionValue(G3D_Map * map, double north, double east, double top,
                    void *value, int type)
{
     int row, col, depth;

     /* convert (north, east, top) into (row, col, depth) */

     row = map->region.rows -
         (north - map->region.south) / (map->region.north -
                                        map->region.south) *
map->region.rows;
     col =
         (east - map->region.west) / (map->region.east -
                                      map->region.west) * map->region.cols;
     depth =
         (top - map->region.bottom) / (map->region.top -
                                       map->region.bottom) *
         map->region.depths;

     /* if (row, col, depth) outside window return NULL value */
     if ((row < 0) || (row >= map->region.rows) ||
         (col < 0) || (col >= map->region.cols) ||
         (depth < 0) || (depth >= map->region.depths)) {
         G3d_setNullValue(value, 1, type);
         return;
     }

     /* get value */
     map->resampleFun(map, row, col, depth, value, type);
}
}}}

r3.out.ascii and r3.in.ascii convert the upper left corner in a lower left
corner, including confusing comments.

r3.in.ascii conversion code:
{{{
for (y = region->rows - 1; y >= 0; y--) /* go south to north */
}}}

r3.out.ascii conversion code:
{{{
for (y = rows - 1; y >= 0; y--) { /* north to south */
}}}

I don't know why it was implemented using this behavior, because this is
not slice compatible to r.out/in.ascii. But the r3.out.ascii manpage says
this:

   One level maps can be imported with r.in.ascii (Raster 2D) after
removing
   the header lines "top", "bottom" and "levels".

which is obviously wrong?

>I also notice that r3.out.vtk exports celldata with the data south-up
within the ascii.vtk file, but I'm not sure if that is part of the VTK
design, or the export module flipping what it doesn't actually have to. ??

This is VTK design. The corner of VTK image data is lower left.

>tangent: I don't know Paraview well enough to reliably say which way
'should' be up there. The first thing I notice in Paraview 3.8.0 after
loading the data, apply, Display|Style->representation,Surface +
Color->Color by->layer name is that the top level is on the bottom and the
bottom layer is on the top. again in the display tab setting
Transformation->Scale->z to -1 makes it ok, but it gets very dark, as if
it is all in shadow. (test data is cube containing 0/1 boolean values,
similar to what r.to.rast3elev would create from a 2.5D DEM for land/air.
I'm still struggling to visualize that well in NVIZ or paraview.) ??or
perhaps I set top: and bottom: backwards in the r3.in.ascii header?
r3.to.rast shows level 0001 as top layer, level 25 (of 25) as bottom
layer. (data is all below sea level)

In case r3.out.vtk is used for Paraview data export, the data should have
a correct coordinate system, no shifting or mirroring is needed.
Maybe the input data is flipped?

Best
Soeren

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/1385#comment:1&gt;
GRASS GIS <http://grass.osgeo.org>

3

#1385: r3.in.ascii man page improvements
-------------------------+--------------------------------------------------
Reporter: hamish | Owner: grass-dev@…
     Type: task | Status: new
Priority: normal | Milestone: 6.4.2
Component: Docs | Version: svn-develbranch6
Keywords: r3.in.ascii | Platform: All
      Cpu: All |
-------------------------+--------------------------------------------------
Changes (by hamish):

  * platform: Linux => All
  * type: defect => task
  * cpu: x86-64 => All

Comment:

Thanks, I was a bit mixed up, both top-bottom, and north-south.

Actually now that I inspect it closely my data was already oriented upside
down with SW corner first (this is Fortran output, maybe that's why the r3
code was that way? ...oops, the fortran code comments even tell me the 'j'
output will inverted. ISTR from matlab weirdness that fortran data arrays
are ordered in memory first by column (not by row), but I'm pretty sure
that is a different matter and irrelevant here)

So yes, r3.in.ascii does as the man page says and wants n-s flipped data,
and the main point of my bug report is invalid-- changing the ticket
description and continuing on as a "improve man page" task.

I'll change the help page example to read:
{{{
,y3,
,y2,
,y1,
,y3,
,y2,
,y1,
...
}}}

and does this sound ok for the r3.in.ascii help page:

{{{
The bottom z-level is given first, the top z-level is given last.
Thus the first data value in the ASCII input file refers to the SW
corner of the bottom-most z-level, and the last data value in the
ASCII file refers to the NE corner of the top-most z-level.
}}}
?

Soeren wrote:
> r3.out.ascii and r3.in.ascii convert the upper left corner
> in a lower left corner
...
> I don't know why it was implemented using this behavior, because
> this is not slice compatible to r.out/in.ascii.

for grass7, should we change that to remove the flip row conversion?
I guess we would add a "version: grass7" to the header for those and issue
an error if the version line is missing? G_fatal_error("Please add a
`version: grass6` line to the header if it is from an earlier version of
GRASS")

> But the r3.out.ascii manpage says this:
> One level maps can be imported with r.in.ascii (Raster 2D) after
> removing the header lines "top", "bottom" and "levels".
>
> which is obviously wrong?

I agree, that must be wrong. (maybe the flipping was added mid-design?)
But we could keep that text in the grass7 version if we decide to
(un)invert the rows there.
[It was a bit hard to understand at first, the word "One" should really be
"Single"]

re. flipped levels:
> In case r3.out.vtk is used for Paraview data export, the data should
> have a correct coordinate system, no shifting or mirroring is needed.
> Maybe the input data is flipped?

my data was giving the top level first, but r3.in.ascii wants the bottom
level first. So r3.out.vtk is fine, PBKAC. I can easily fix that with
r3.to.rast and then r.to.rast3 with the level map order reversed.

----
some more questions / ideas for grass7 before cleaning up the r3.in.ascii
man page:

  * re. the type={float|double} option -- there is no CELL version for G3D
rasters?

  * precision: default, max, 0-52
  -- if double is used, any point to go beyond %.15g (or %.16g) ?

  for *._in_.ascii why does it matter? shouldn't that be a function of
type=CELL,FCELL,DCELL? shouldn't that option just be applicable to
*.out.ascii?

  * nv: make the default "*" to match other grass ascii module convention?

  * compression: I can't remember, and r.compress doesn't mention it:
   was lzw ever reenabled after the patent expired?

  * tiledimension: this does override the values in the file header? does
it obviate the need to give header values for rows,cols,levels?

cheers,
Hamish

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/1385#comment:2&gt;
GRASS GIS <http://grass.osgeo.org>

4

#1385: r3.in.ascii man page improvements
-------------------------+--------------------------------------------------
Reporter: hamish | Owner: grass-dev@…
     Type: task | Status: new
Priority: normal | Milestone: 6.4.2
Component: Docs | Version: svn-develbranch6
Keywords: r3.in.ascii | Platform: All
      Cpu: All |
-------------------------+--------------------------------------------------

Comment(by huhabla):

Replying to [comment:2 hamish]:
...
> I'll change the help page example to read:
> {{{
> ,y3,
> ,y2,
> ,y1,
> ,y3,
> ,y2,
> ,y1,
> ...
> }}}
>
> and does this sound ok for the r3.in.ascii help page:
>
> {{{
> The bottom z-level is given first, the top z-level is given last.
> Thus the first data value in the ASCII input file refers to the SW
> corner of the bottom-most z-level, and the last data value in the
> ASCII file refers to the NE corner of the top-most z-level.
> }}}
> ?

Thats fine.

>
> Soeren wrote:
> > r3.out.ascii and r3.in.ascii convert the upper left corner
> > in a lower left corner
> ...
> > I don't know why it was implemented using this behavior, because
> > this is not slice compatible to r.out/in.ascii.
>
> for grass7, should we change that to remove the flip row conversion?
> I guess we would add a "version: grass7" to the header for those and
issue an error if the version line is missing? G_fatal_error("Please add
a `version: grass6` line to the header if it is from an earlier version of
GRASS")

I think we should support several different formats of data ordering in
r3.out.ascii and r3.in.ascii:
   * north-south bottom-top order as default (nsbt)
   * south-north bottom-top order (snbt) -> the current standard
   * south-north top-bottom order (sntb)

A new header entry "order:nsbt" will specify the storage structure. So two
new leading header lines can be added:

{{{
version:grass7
order:nsbt
north:...
south:...
}}}
In case the new header lines are not present, the old order scheme is
assumed and imported by r3.in.ascii. This will make it backward
compatible. A new option can be added to r3.out.ascii: "order" with values
"nsbt, snbt, sntb" and the flag -b to create "backward compatible" grass6
ascii output using the snbt order without the new header lines.

What do you think?

>
>
> > But the r3.out.ascii manpage says this:
> > One level maps can be imported with r.in.ascii (Raster 2D) after
> > removing the header lines "top", "bottom" and "levels".
> >
> > which is obviously wrong?
>
> I agree, that must be wrong. (maybe the flipping was added mid-design?)
But we could keep that text in the grass7 version if we decide to
(un)invert the rows there.
> [It was a bit hard to understand at first, the word "One" should really
be "Single"]
>
>
> re. flipped levels:
> > In case r3.out.vtk is used for Paraview data export, the data should
> > have a correct coordinate system, no shifting or mirroring is needed.
> > Maybe the input data is flipped?
>
> my data was giving the top level first, but r3.in.ascii wants the bottom
level first. So r3.out.vtk is fine, PBKAC. I can easily fix that with
r3.to.rast and then r.to.rast3 with the level map order reversed.
>
>
> ----
> some more questions / ideas for grass7 before cleaning up the
r3.in.ascii man page:
>
>
> * re. the type={float|double} option -- there is no CELL version for
G3D rasters?#

No, the type CELL is not present. Only FCELL and DCELL.

>
>
> * precision: default, max, 0-52
> -- if double is used, any point to go beyond %.15g (or %.16g) ?
>
> for *._in_.ascii why does it matter? shouldn't that be a function of
type=CELL,FCELL,DCELL? shouldn't that option just be applicable to
*.out.ascii?

The precision value is the number of digits used as mantissa in the
internal storage 0 -23 for float and 0 - 52 for double -> this is only
relevant for compression and NULL data handling.

>
>
> * nv: make the default "*" to match other grass ascii module
convention?

Good idea.

>
>
> * compression: I can't remember, and r.compress doesn't mention it:
> was lzw ever reenabled after the patent expired?

No clue.

>
>
> * tiledimension: this does override the values in the file header? does
it obviate the need to give header values for rows,cols,levels?

No, the tile size is used for internal storage and does not affect any
header information. G3d uses a (cached) tile based approach rather then a
row based approach. This allows arbitrary value access when tile caching
is enabled ... which is really cool.

Best
Soeren

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/1385#comment:3&gt;
GRASS GIS <http://grass.osgeo.org>

5

#1385: r3.in.ascii man page improvements
-------------------------+--------------------------------------------------
Reporter: hamish | Owner: grass-dev@…
     Type: task | Status: new
Priority: normal | Milestone: 6.4.2
Component: Docs | Version: svn-develbranch6
Keywords: r3.in.ascii | Platform: All
      Cpu: All |
-------------------------+--------------------------------------------------

Comment(by hamish):

Hi,

did you find the r3.to.rast -> 0 values for FCELL rast3d maps bug already?
or should I open a new report for that?

re. output precision, anything more than %.16g or %.17g for IEEE double
precision FPs is just noise. There's no point of going to 52.

any documented reason why there is no CELL support for G3D?

for r.out.ascii I would like to set the default output precision
dynamically (%.7g for FCELL, %.15g for DCELL, [%ld for CELL??]). Maybe we
could try that in gr7 and see how it looks?
(I just added a quick&rough r.flip addon script using that, and the round
trip is lossy; Yann's upcoming? C port of grass5's version won't have that
problem)

thanks,
Hamish

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/1385#comment:4&gt;
GRASS GIS <http://grass.osgeo.org>

6

#1385: r3.in.ascii man page improvements
-------------------------+--------------------------------------------------
Reporter: hamish | Owner: grass-dev@…
     Type: task | Status: new
Priority: normal | Milestone: 6.4.2
Component: Docs | Version: svn-develbranch6
Keywords: r3.in.ascii | Platform: All
      Cpu: All |
-------------------------+--------------------------------------------------

Comment(by hamish):

> The precision value is the number of digits used as mantissa in the
> internal storage 0 -23 for float and 0 - 52 for double -> this is only
> relevant for compression and NULL data handling.

ah, ok.

(I think we should expand tile and precision option descriptions+labels to
clarify this)

Hamish

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/1385#comment:5&gt;
GRASS GIS <http://grass.osgeo.org>

7

#1385: r3.in.ascii man page improvements
-------------------------+--------------------------------------------------
Reporter: hamish | Owner: grass-dev@…
     Type: task | Status: new
Priority: normal | Milestone: 6.4.2
Component: Docs | Version: svn-develbranch6
Keywords: r3.in.ascii | Platform: All
      Cpu: All |
-------------------------+--------------------------------------------------

Comment(by huhabla):

Replying to [comment:4 hamish]:
> Hi,
>
> did you find the r3.to.rast -> 0 values for FCELL rast3d maps bug
already? or should I open a new report for that?

Uh? I did not know that such a bug exists. Can you please open a ticket
and provide a synthetic data test?

> any documented reason why there is no CELL support for G3D?

Unfortunately not. I did not found any information why CELL value support
was not implemented. I guess Helena may know why?

>
> for r.out.ascii I would like to set the default output precision
dynamically (%.7g for FCELL, %.15g for DCELL, [%ld for CELL??]). Maybe we
could try that in gr7 and see how it looks?
> (I just added a quick&rough r.flip addon script using that, and the
round trip is lossy; Yann's upcoming? C port of grass5's version won't
have that problem)

Feel free to modify the code. :slight_smile:

Best regards
Soeren

>
>
> thanks,
> Hamish

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/1385#comment:6&gt;
GRASS GIS <http://grass.osgeo.org>

8

#1385: r3.in.ascii man page improvements
-------------------------+--------------------------------------------------
Reporter: hamish | Owner: grass-dev@…
     Type: task | Status: new
Priority: normal | Milestone: 6.4.2
Component: Raster3D | Version: svn-develbranch6
Keywords: r3.in.ascii | Platform: All
      Cpu: All |
-------------------------+--------------------------------------------------
Changes (by hamish):

  * component: Docs => Raster3D

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/1385#comment:7&gt;
GRASS GIS <http://grass.osgeo.org>