[GRASS5] [bug #992] (grass) g.region Vask problem in lat/long

Justin_Hickey_via_RT · March 27, 2002, 5:10pm

this bug's URL: http://intevation.de/rt/webrt?serial_num=992
-------------------------------------------------------------------------

Subject: g.region Vask problem in lat/long

Platform: GNU/Linux/i386
grass obtained from: CVS
grass binary for platform: Compiled from Sources

Working in lat/long location, I found a small (unusual) error:

GRASS:~ > g.region -p
projection: 3 (Latitude-Longitude)
zone: 0
datum: WGS84
ellipsoid: WGS84
north: 17:46S
south: 26:54:30S
west: 19:59E
east: 29:22E
nsres: 0:00:30
ewres: 0:00:30
rows: 1097
cols: 1126

Now starting g.region, selecting menu 1:
V_ask error: Length out of bounds in call to V_const <- this warning appears for 3 sec

IDENTIFY REGION

                   Default GRID RESOLUTION Region
                     0:00:30--- East-West --- 0:00:30___
                     0:00:30-- North-South -- 0:00:30___

The west value makes me unhappy. Looks like a rounding problem somewhere.

Markus

-------------------------------------------- Managed by Request Tracker

Glynn_Clements · March 27, 2002, 6:26pm

Request Tracker wrote:

Subject: g.region Vask problem in lat/long

Platform: GNU/Linux/i386
grass obtained from: CVS
grass binary for platform: Compiled from Sources

Working in lat/long location, I found a small (unusual) error:

GRASS:~ > g.region -p

Can you tell us the default region (rather than the current region),
as it seems to be the default West value that's causing the problem.

Now starting g.region, selecting menu 1:
V_ask error: Length out of bounds in call to V_const <- this warning appears for 3 sec

The west value makes me unhappy. Looks like a rounding problem somewhere.

Well, floating point can't represent exactly fractions where the
denominator isn't a power of two, so rounding errors will occur.

G_lon_format() (which generates the string in question) doesn't use
any specific precision, so you get the standard behaviour of "%f",
which is 6 digits after the point.

But how many digits should it use? Are whole seconds sufficiently
accurate? Unless you require accuracy better than thousandths of a
second, the argument should be converted to seconds and integer
arithmetic used.

--
Glynn Clements <glynn.clements@virgin.net>

Markus_Neteler1 · March 27, 2002, 9:50pm

On Wed, Mar 27, 2002 at 06:26:07PM +0000, Glynn Clements wrote:

Request Tracker wrote:

> Subject: g.region Vask problem in lat/long
>
> Platform: GNU/Linux/i386
> grass obtained from: CVS
> grass binary for platform: Compiled from Sources
>
> Working in lat/long location, I found a small (unusual) error:
>
> GRASS:~ > g.region -p

Can you tell us the default region (rather than the current region),
as it seems to be the default West value that's causing the problem.

Yes, please have a look in the RT - I had submitted a followup (and
obviously forgot to cc to the list). In fact the default West value that's
causing the problem.

> Now starting g.region, selecting menu 1:
> V_ask error: Length out of bounds in call to V_const <- this warning appears for 3 sec

> The west value makes me unhappy. Looks like a rounding problem somewhere.

Well, floating point can't represent exactly fractions where the
denominator isn't a power of two, so rounding errors will occur.

Mhhh, but I was demanding rounded coordinates with
g.region res=0:00:30 -a
earlier. The DEFAULT_WIND contains a round West, only inside g.region
it appears strange (means: should be everywhere the same value).

G_lon_format() (which generates the string in question) doesn't use
any specific precision, so you get the standard behaviour of "%f",
which is 6 digits after the point.

But how many digits should it use? Are whole seconds sufficiently
accurate? Unless you require accuracy better than thousandths of a
second, the argument should be converted to seconds and integer
arithmetic used.

Eventually yes. At least the overwriting (see RT followup) isn't nice.

Markus

Glynn_Clements · March 27, 2002, 10:50pm

Markus Neteler wrote:

> Can you tell us the default region (rather than the current region),
> as it seems to be the default West value that's causing the problem.

Yes, please have a look in the RT - I had submitted a followup (and
obviously forgot to cc to the list). In fact the default West value that's
causing the problem.

OK, I've looked at that, but I cannot reproduce the problem. Can you
send me the files?

> > Now starting g.region, selecting menu 1:
> > V_ask error: Length out of bounds in call to V_const <- this warning appears for 3 sec
>
> > The west value makes me unhappy. Looks like a rounding problem somewhere.
>
> Well, floating point can't represent exactly fractions where the
> denominator isn't a power of two, so rounding errors will occur.

Mhhh, but I was demanding rounded coordinates with
g.region res=0:00:30 -a
earlier. The DEFAULT_WIND contains a round West, only inside g.region
it appears strange (means: should be everywhere the same value).

Well, GRASS doesn't use DMS internally; the DMS -> FP -> DMS
conversions could easily introduce rounding error. Particularly if
something is multiplying a rounded representation of one second by 30.

> G_lon_format() (which generates the string in question) doesn't use
> any specific precision, so you get the standard behaviour of "%f",
> which is 6 digits after the point.
>
> But how many digits should it use? Are whole seconds sufficiently
> accurate? Unless you require accuracy better than thousandths of a
> second, the argument should be converted to seconds and integer
> arithmetic used.

Eventually yes. At least the overwriting (see RT followup) isn't nice.

    I slighly change the report:
    The field of default region in the g.region menu is too long and writing
    into the current region. If there are no objections we should limit the
    field length in
    src/libes/gis/edit_cellhd.c

What should it be changed to? To whatever space is currently
available?

This would be better done in G_{lat,lon}_format(); using a precision
specifier (e.g. "%.2f") will round the result, so that e.g.
"59.999999" -> "60.00", while truncating the string within
G_edit_cellhd() would, er, truncate it, so "59.999999" -> "59.99".

But G_{lat,lon}_format() are used elsewhere, and other uses might
warrant greater precision. Ideally the formatting functions should
allow the precision to be specified, but that has compatibility
implications.

Ultimately, there are two problems here. One is that the error is
being introduced in the first place, the other is that G_edit_cellhd()
only copes with certain easting/northing values.

--
Glynn Clements <glynn.clements@virgin.net>

Markus_Neteler1 · March 28, 2002, 7:59am

On Wed, Mar 27, 2002 at 10:50:50PM +0000, Glynn Clements wrote:

> > Markus Neteler wrote:
> > > Now starting g.region, selecting menu 1:
> > > V_ask error: Length out of bounds in call to V_const <- this warning appears for 3 sec
> >
> > > The west value makes me unhappy. Looks like a rounding problem somewhere.
> >
> > Well, floating point can't represent exactly fractions where the
> > denominator isn't a power of two, so rounding errors will occur.
>
> Mhhh, but I was demanding rounded coordinates with
> g.region res=0:00:30 -a
> earlier. The DEFAULT_WIND contains a round West, only inside g.region
> it appears strange (means: should be everywhere the same value).

I am sorry, Glynn - today I could only partially reproduce the problem.
It seems I was overlooking something yesterday.

g.region -dp
projection: 3 (Latitude-Longitude)
zone: 0
datum: WGS84
ellipsoid: WGS84
north: 8:59:45S
south: 29:00:15S
west: 15:59:59.999982E
east: 34:00:29.999982E
nsres: 0:00:30
ewres: 0:00:30
rows: 2401
cols: 2161

Now the error
"V_ask error: Length out of bounds in call to V_const"
is explanable.

I generated this location from GLOBE DEM30 (user defined window)
automatically with r.in.gdal (it reads GLOBE DEM30 directly),
therefore these strange west/east boundaries:

gdalinfo botswana.bil
Driver: EHdr/ESRI .hdr Labelled
Size is 2160, 1920
Coordinate System is '
Origin = (16.000000,-13.000000)
Pixel Size = (0.008333,-0.008333)
Corner Coordinates:
Upper Left ( 16.0000000, -13.0000000)
Lower Left ( 16.0000000, -28.9999936)
Upper Right ( 33.9999928, -13.0000000)
Lower Right ( 33.9999928, -28.9999936)
Center ( 24.9999964, -20.9999968)
Band 1 Block=2160x1 Type=Int16, ColorInterp=Undefined
NoData Value=-500

Well, GRASS doesn't use DMS internally; the DMS -> FP -> DMS
conversions could easily introduce rounding error. Particularly if
something is multiplying a rounded representation of one second by 30.

> > G_lon_format() (which generates the string in question) doesn't use
> > any specific precision, so you get the standard behaviour of "%f",
> > which is 6 digits after the point.
> >
> > But how many digits should it use? Are whole seconds sufficiently
> > accurate? Unless you require accuracy better than thousandths of a
> > second, the argument should be converted to seconds and integer
> > arithmetic used.
>
> Eventually yes. At least the overwriting (see RT followup) isn't nice.

    I slighly change the report:
    The field of default region in the g.region menu is too long and writing
    into the current region. If there are no objections we should limit the
    field length in
    src/libes/gis/edit_cellhd.c

What should it be changed to? To whatever space is currently
available?

I think so. Or we have to use two lines

       |
default_west
       | current_west
       |
In this case default_west may keep its original length. I am not
sure if that is a good idea.

This would be better done in G_{lat,lon}_format(); using a precision
specifier (e.g. "%.2f") will round the result, so that e.g.
"59.999999" -> "60.00", while truncating the string within
G_edit_cellhd() would, er, truncate it, so "59.999999" -> "59.99".

But G_{lat,lon}_format() are used elsewhere, and other uses might
warrant greater precision. Ideally the formatting functions should
allow the precision to be specified, but that has compatibility
implications.

Ultimately, there are two problems here. One is that the error is
being introduced in the first place, the other is that G_edit_cellhd()
only copes with certain easting/northing values.

If still needed, I can put the entire location somewhere. Let me know
(46MB uncompressed).

Markus

Glynn_Clements · March 28, 2002, 12:04pm

Markus Neteler wrote:

I generated this location from GLOBE DEM30 (user defined window)
automatically with r.in.gdal (it reads GLOBE DEM30 directly),
therefore these strange west/east boundaries:

gdalinfo botswana.bil
Driver: EHdr/ESRI .hdr Labelled
Size is 2160, 1920
Coordinate System is '
Origin = (16.000000,-13.000000)
Pixel Size = (0.008333,-0.008333)
Corner Coordinates:
Upper Left ( 16.0000000, -13.0000000)
Lower Left ( 16.0000000, -28.9999936)
Upper Right ( 33.9999928, -13.0000000)
Lower Right ( 33.9999928, -28.9999936)
Center ( 24.9999964, -20.9999968)
Band 1 Block=2160x1 Type=Int16, ColorInterp=Undefined
NoData Value=-500

Well, that explains where the error is coming from. Perhaps you should
correspond with Frank on this aspect.

> > Eventually yes. At least the overwriting (see RT followup) isn't nice.
>
> I slighly change the report:
> The field of default region in the g.region menu is too long and writing
> into the current region. If there are no objections we should limit the
> field length in
> src/libes/gis/edit_cellhd.c
>
> What should it be changed to? To whatever space is currently
> available?

I think so. Or we have to use two lines

       |
default_west
       | current_west
       |
In this case default_west may keep its original length. I am not
sure if that is a good idea.

I'll look into rounding the data before passing it to the formatting
routines. I'll also truncate the string to prevent errors.

If still needed, I can put the entire location somewhere. Let me know
(46MB uncompressed).

Not necessary.

--
Glynn Clements <glynn.clements@virgin.net>

Markus_Neteler1 · March 28, 2002, 1:41pm

On Thu, Mar 28, 2002 at 12:04:17PM +0000, Glynn Clements wrote:

Markus Neteler wrote:

> I generated this location from GLOBE DEM30 (user defined window)
> automatically with r.in.gdal (it reads GLOBE DEM30 directly),
> therefore these strange west/east boundaries:
>
> gdalinfo botswana.bil
> Driver: EHdr/ESRI .hdr Labelled
> Size is 2160, 1920
> Coordinate System is '
> Origin = (16.000000,-13.000000)
> Pixel Size = (0.008333,-0.008333)
> Corner Coordinates:
> Upper Left ( 16.0000000, -13.0000000)
> Lower Left ( 16.0000000, -28.9999936)
> Upper Right ( 33.9999928, -13.0000000)
> Lower Right ( 33.9999928, -28.9999936)
> Center ( 24.9999964, -20.9999968)
> Band 1 Block=2160x1 Type=Int16, ColorInterp=Undefined
> NoData Value=-500

Well, that explains where the error is coming from. Perhaps you should
correspond with Frank on this aspect.

To reduce confusion:
Digging in my mind I remember that I was creating this location
eventually from a different source. The second map (GLOBE 30) then
extended the location. So the g.region and gdalinfo must not necessarily
match.

To my experience the r.in.gdal runs well.

Thanks for looking into g.region,

Markus

Glynn_Clements · March 28, 2002, 7:23pm

Glynn Clements wrote:

I'll look into rounding the data before passing it to the formatting
routines. I'll also truncate the string to prevent errors.

Done.

The values displayed by G_edit_cellhd() are rounded to fit into 10
characters; also the result is truncated, just in case the rounding
doesn't do the job.

This provides a resolution of 1 second for lat-lon, and 1 unit (metre,
foot) for other projections. It might be possible to allow greater
resolution for non-lat/lon projections; it depends upon the maximum
range of the values.

--
Glynn Clements <glynn.clements@virgin.net>