r3.in.ascii man page improvements

[hamish]

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 :)

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.

[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

[hamish]

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

[huhabla]

Replying to 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

[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

[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

[huhabla]

Replying to 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. :)

Best regards Soeren

thanks, Hamish

[neteler]

Anything left open with respect to the original report?