[GRASS5] r.in.gdal - precision problem in lib/gis/adj_cellhd.c

GRASS GIS GRASS GIS developer
,
1

Hi,

I tried to import the Europe tile of SRTM/GTOPO30 from GLCF.
There is a small issue which makes GRASS refuse the file:

gdalinfo SRTM_GTOPO_u30_n090w020.tif
Driver: GTiff/GeoTIFF
Size is 4800, 6000
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.2572235629972,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["metre",1],
    AUTHORITY["EPSG","4326"]]
Origin = (-20.000000,90.000000)
Pixel Size = (0.00833333,-0.00833333)
Metadata:
  TIFFTAG_DOCUMENTNAME=D:\Kuan\srtm\SRTM30\w020n90_new.tif
  TIFFTAG_IMAGEDESCRIPTION=IDL TIFF file
  TIFFTAG_XRESOLUTION=100
  TIFFTAG_YRESOLUTION=100
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
Corner Coordinates:
Upper Left ( -20.0000000, 90.0000000) ( 20d 0'0.00"W, 90d 0'0.00"N)
Lower Left ( -20.0000000, 39.9999974) ( 20d 0'0.00"W, 39d59'59.99"N)
Upper Right ( 20.0000021, 90.0000000) ( 20d 0'0.01"E, 90d 0'0.00"N)
Lower Right ( 20.0000021, 39.9999974) ( 20d 0'0.01"E, 39d59'59.99"N)
Center ( 0.0000010, 64.9999987) ( 0d 0'0.00"E, 65d 0'0.00"N)
Band 1 Block=4800x1 Type=Float32, ColorInterp=Gray

-> looks reasonable

GRASS 6.1.cvs (latlong):~ > r.in.gdal SRTM_GTOPO_u30_n090w020.tif out=srtm90_gtopo30
A datum name wgs84 (WGS_1984) was specified without transformation parameters.
Note that the GRASS default for wgs84 is towgs84=0.000,0.000,0.000.
Projection of input dataset and current location appear to match.
Proceeding with import...
D0/0: North: 90.000000
WARNING: G_set_window(): Illegal latitude for North

-> I have added debug output to lib/gis/adj_cellhd.c
--> strange, but (using %.10f output):

GRASS 6.1.cvs (latlong):~ > r.in.gdal SRTM_GTOPO_u30_n090w020.tif out=srtm90_gtopo30
A datum name wgs84 (WGS_1984) was specified without transformation parameters.
Note that the GRASS default for wgs84 is towgs84=0.000,0.000,0.000.
Projection of input dataset and current location appear to match.
Proceeding with import...
D0/0: North: 90.0000000002
WARNING: G_set_window(): Illegal latitude for North

-> to me it looks like a C problem.

How can I convince lib/gis/adj_cellhd.c to work as desired?
The comparison is currently:

if (cellhd->north > 90.0)
     return (_("Illegal latitude for North"));

Thanks,
Markus

2

As a followup:

Hacking lib/gis/adj_cellhd.c to permit the import,
I get the following cell header:

grassdata/latlong/PERMANENT/cellhd/srtm90_gtopo30
proj: 3
zone: 0
north: 90:00:00.000001N
south: 39:59:59.990613N
east: 20:00:00.007509E
west: 20:00:00.000001W
cols: 4800
rows: 6000
e-w resol: 0:00:30.000002
n-s resol: 0:00:30.000002
format: -1
compressed: 2

Probably it's a GRASS/GDAL issue with LatLong data?

Fixing the cell header files manually solved the
problem but isn't very elegant.

Any hints are welcome

Markus

On Wed, Apr 06, 2005 at 03:23:00PM +0200, Markus Neteler wrote:

Hi,

I tried to import the Europe tile of SRTM/GTOPO30 from GLCF.
There is a small issue which makes GRASS refuse the file:

gdalinfo SRTM_GTOPO_u30_n090w020.tif
Driver: GTiff/GeoTIFF
Size is 4800, 6000
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.2572235629972,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["metre",1],
    AUTHORITY["EPSG","4326"]]
Origin = (-20.000000,90.000000)
Pixel Size = (0.00833333,-0.00833333)
Metadata:
  TIFFTAG_DOCUMENTNAME=D:\Kuan\srtm\SRTM30\w020n90_new.tif
  TIFFTAG_IMAGEDESCRIPTION=IDL TIFF file
  TIFFTAG_XRESOLUTION=100
  TIFFTAG_YRESOLUTION=100
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
Corner Coordinates:
Upper Left ( -20.0000000, 90.0000000) ( 20d 0'0.00"W, 90d 0'0.00"N)
Lower Left ( -20.0000000, 39.9999974) ( 20d 0'0.00"W, 39d59'59.99"N)
Upper Right ( 20.0000021, 90.0000000) ( 20d 0'0.01"E, 90d 0'0.00"N)
Lower Right ( 20.0000021, 39.9999974) ( 20d 0'0.01"E, 39d59'59.99"N)
Center ( 0.0000010, 64.9999987) ( 0d 0'0.00"E, 65d 0'0.00"N)
Band 1 Block=4800x1 Type=Float32, ColorInterp=Gray

-> looks reasonable

GRASS 6.1.cvs (latlong):~ > r.in.gdal SRTM_GTOPO_u30_n090w020.tif out=srtm90_gtopo30
A datum name wgs84 (WGS_1984) was specified without transformation parameters.
Note that the GRASS default for wgs84 is towgs84=0.000,0.000,0.000.
Projection of input dataset and current location appear to match.
Proceeding with import...
D0/0: North: 90.000000
WARNING: G_set_window(): Illegal latitude for North

-> I have added debug output to lib/gis/adj_cellhd.c
--> strange, but (using %.10f output):

GRASS 6.1.cvs (latlong):~ > r.in.gdal SRTM_GTOPO_u30_n090w020.tif out=srtm90_gtopo30
A datum name wgs84 (WGS_1984) was specified without transformation parameters.
Note that the GRASS default for wgs84 is towgs84=0.000,0.000,0.000.
Projection of input dataset and current location appear to match.
Proceeding with import...
D0/0: North: 90.0000000002
WARNING: G_set_window(): Illegal latitude for North

-> to me it looks like a C problem.

How can I convince lib/gis/adj_cellhd.c to work as desired?
The comparison is currently:

if (cellhd->north > 90.0)
     return (_("Illegal latitude for North"));

Thanks,
Markus

3

On Wed, 2005-04-06 at 15:23 +0200, Markus Neteler wrote:

How can I convince lib/gis/adj_cellhd.c to work as desired?
The comparison is currently:

if (cellhd->north > 90.0)
     return (_("Illegal latitude for North"));

float north;

sscanf(cellhd->north, "%.1f", &north);
if (north > 90.0)
    return (_("Illegal latitude for North.\n"));

--
Brad Douglas <rez@touchofmadness.com>

4

Hacking lib/gis/adj_cellhd.c to permit the import,
I get the following cell header:

grassdata/latlong/PERMANENT/cellhd/srtm90_gtopo30
proj: 3
zone: 0
north: 90:00:00.000001N
south: 39:59:59.990613N
east: 20:00:00.007509E
west: 20:00:00.000001W
cols: 4800
rows: 6000
e-w resol: 0:00:30.000002
n-s resol: 0:00:30.000002
format: -1
compressed: 2

Probably it's a GRASS/GDAL issue with LatLong data?

Fixing the cell header files manually solved the
problem but isn't very elegant.

Any hints are welcome

It's not trying to add half the cell to the top row of the TIFF, thus
going beyond 90d 0' 15"?

I thought SRTM only made it to 60 deg lat anyway? Space Shuttle is not
polar orbiting.

I have noticed in the past that v.out.ogr to a shapefile does similar
strange off-by 1e-9 type errors to double precision columns in the
data.

Hamish

5

Brad Douglas wrote:

On Wed, 2005-04-06 at 15:23 +0200, Markus Neteler wrote:

How can I convince lib/gis/adj_cellhd.c to work as desired?
The comparison is currently:

if (cellhd->north > 90.0)
    return (_("Illegal latitude for North"));

float north;

sscanf(cellhd->north, "%.1f", &north);
if (north > 90.0)
    return (_("Illegal latitude for North.\n"));

It can silently damage data. I am not in favour.

Radim

6

On Thu, Apr 07, 2005 at 04:25:38PM +1200, Hamish wrote:

> Hacking lib/gis/adj_cellhd.c to permit the import,
> I get the following cell header:
>
> grassdata/latlong/PERMANENT/cellhd/srtm90_gtopo30
> proj: 3
> zone: 0
> north: 90:00:00.000001N
> south: 39:59:59.990613N
> east: 20:00:00.007509E
> west: 20:00:00.000001W
> cols: 4800
> rows: 6000
> e-w resol: 0:00:30.000002
> n-s resol: 0:00:30.000002
> format: -1
> compressed: 2
>
> Probably it's a GRASS/GDAL issue with LatLong data?
>
> Fixing the cell header files manually solved the
> problem but isn't very elegant.
>
> Any hints are welcome

It's not trying to add half the cell to the top row of the TIFF, thus
going beyond 90d 0' 15"?

Probably. But then we would exclude a lot of global data!

I thought SRTM only made it to 60 deg lat anyway? Space Shuttle is not
polar orbiting.

Note that these data are the GTOPO30 data enhanced by SRTM, not
the original SRTM data. So north of 60 deg GTOPO30 remained untouched AFAIK
(northern hemisphere).

I have noticed in the past that v.out.ogr to a shapefile does similar
strange off-by 1e-9 type errors to double precision columns in the
data.

Mhh, a systematic problem?

Markus

7

On Thu, 2005-04-07 at 08:59 +0200, Radim Blazek wrote:

Brad Douglas wrote:
> On Wed, 2005-04-06 at 15:23 +0200, Markus Neteler wrote:
>
>>How can I convince lib/gis/adj_cellhd.c to work as desired?
>>The comparison is currently:
>>
>> if (cellhd->north > 90.0)
>> return (_("Illegal latitude for North"));
>
>
> float north;
>
> sscanf(cellhd->north, "%.1f", &north);
> if (north > 90.0)
> return (_("Illegal latitude for North.\n"));
>

It can silently damage data. I am not in favour.

It should suffice for this specific purpose as a short-term fix.

Although shunned upon, that usage of sscanf() is littered throughout
GRASS code. I say we temporarily commit the above to get it "mostly"
working and work on a permanent solution to non-integer comparisons.

Maybe we could do one of these:

/* Compares two floats to 'precision' decimal places */
/* Returns 1 if equal, 0 otherwise */
int G_f_equal(float a, float b, int precision)
...

or:

/* Compares two floats; Precision "adjusted" */
/* Returns 1 if equal, 0 otherwise */
int G_relequal(float a, float b)
{
#define ABS(x) ((x) < 0 ? -(x) : (x))
#define MAX(x, y) ((x) <= (y) ? (y) : (x))

    float c = ABS(a);
    float d = ABS(b);

    c = MAX(c, d);
    d = ABS(a - b) / c;

   return (d == 0. ? 1 : 0);
}

/* Compares two floats. "Precision adjusted" */
/* Returns positive, 0, or negative value if a */
/* is greater, equal, or less than b */
float G_relcompare(float a, float b)
...

Do any of these sound useful?

--
Brad Douglas <rez@touchofmadness.com>

8

On Apr 7, 2005 11:09 AM, Markus Neteler <neteler@itc.it> wrote:

> It's not trying to add half the cell to the top row of the TIFF, thus
> going beyond 90d 0' 15"?

Probably. But then we would exclude a lot of global data!

Markus,

It seems like a precision issue but it is hard to know for use since stuff
like gdalinfo does not print out the origin with full double precision values.
I am not positive whether the precision issue is in the original GeoTIFF
file, GDAL, v.in.gdal or what. I suspect the original file.

It might pay to have adj_cellhd() to some subtle fixing up in cases that
just exceed global bounds by some small epsilon.

Best regards,
--
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up | Frank Warmerdam, warmerdam@pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush | Geospatial Programmer for Rent

9

> I have noticed in the past that v.out.ogr to a shapefile does
> similar strange off-by 1e-9 type errors to double precision columns
> in the data.

Mhh, a systematic problem?

No I think it is probably separate but wanted to note. A double cast to
float and then back to double introducing some noise?

Hamish

10

Brad Douglas wrote:

sscanf(cellhd->north, "%.1f", &north);
if (north > 90.0)
   return (_("Illegal latitude for North.\n"));

It can silently damage data. I am not in favour.

It should suffice for this specific purpose as a short-term fix.

It does not fix the problem, it just hides it. If we do that, nobody will fix that in next 10 years.

Although shunned upon, that usage of sscanf() is littered throughout
GRASS code. I say we temporarily commit the above to get it "mostly"
working and work on a permanent solution to non-integer comparisons.

Maybe we could do one of these:

/* Compares two floats to 'precision' decimal places */
/* Returns 1 if equal, 0 otherwise */
int G_f_equal(float a, float b, int precision)

/* Compares two floats. "Precision adjusted" */
/* Returns positive, 0, or negative value if a */
/* is greater, equal, or less than b */
float G_relcompare(float a, float b)

But the result will be a map with north > 90 in any case.

If we need a short term solution, I suggest to change the error to warning or maybe a warning for 90 < north < 90.000x and error for north > 90.000x?

Radim

11

But the result will be a map with north > 90 in any case.

If we need a short term solution, I suggest to change the error to
warning or maybe a warning for 90 < north < 90.000x and error for
north
> 90.000x?

Can we establish from that 30" data causing the problem if the top cells
are centered at 90N or at 89d 59' 45" ? Having a row of (non-uniform)
data centered at 90N seems like a mistake in the data creation to me.

I know the ETOPO2 dataset correctly has top edge at 90N after import,
center of top cells are at 89d 59".

Hamish

12

On Sun, Apr 10, 2005 at 02:44:17PM +1200, Hamish wrote:

> But the result will be a map with north > 90 in any case.
>
> If we need a short term solution, I suggest to change the error to
> warning or maybe a warning for 90 < north < 90.000x and error for
> north
> > 90.000x?

Can we establish from that 30" data causing the problem if the top cells
are centered at 90N or at 89d 59' 45" ? Having a row of (non-uniform)
data centered at 90N seems like a mistake in the data creation to me.

As far as I know the GDAL software reports cell *corners* including
the C/C++ functions:

gdalinfo SRTM_GTOPO_u30_n090w020.tif
Driver: GTiff/GeoTIFF
Size is 4800, 6000
...
Origin = (-20.000000,90.000000)
Pixel Size = (0.00833333,-0.00833333)
...
Corner Coordinates:
Upper Left ( -20.0000000, 90.0000000) ( 20d 0'0.00"W, 90d 0'0.00"N)
Lower Left ( -20.0000000, 39.9999974) ( 20d 0'0.00"W, 39d59'59.99"N)
Upper Right ( 20.0000021, 90.0000000) ( 20d 0'0.01"E, 90d 0'0.00"N)
Lower Right ( 20.0000021, 39.9999974) ( 20d 0'0.01"E, 39d59'59.99"N)
Center ( 0.0000010, 64.9999987) ( 0d 0'0.00"E, 65d 0'0.00"N)

The number of rows/columns looks ok.
So GRASS north coordinates should be GDAL corners minus 0.5 resolution.
Or not?

Even if the data set was not generated properly, GRASS should somehow
deal with that (as I cannot ask GLCF to redo the file).

What about a subtle fixing up in cases that just the coordinates
exceed the global bounds by some small epsilon?

Markus

13

(cc grass5)

On Mon, Apr 11, 2005 at 03:12:44PM +0200, Morten Hulden wrote:

Markus Neteler wrote:

>gdalinfo SRTM_GTOPO_u30_n090w020.tif
>Driver: GTiff/GeoTIFF
>Size is 4800, 6000
>...
>Origin = (-20.000000,90.000000)
>Pixel Size = (0.00833333,-0.00833333)
>...
>Corner Coordinates:
>Upper Left ( -20.0000000, 90.0000000) ( 20d 0'0.00"W, 90d 0'0.00"N)
>Lower Left ( -20.0000000, 39.9999974) ( 20d 0'0.00"W, 39d59'59.99"N)
>Upper Right ( 20.0000021, 90.0000000) ( 20d 0'0.01"E, 90d 0'0.00"N)
>Lower Right ( 20.0000021, 39.9999974) ( 20d 0'0.01"E, 39d59'59.99"N)
>Center ( 0.0000010, 64.9999987) ( 0d 0'0.00"E, 65d 0'0.00"N)
>
>The number of rows/columns looks ok.
>So GRASS north coordinates should be GDAL corners minus 0.5 resolution.
>Or not?
>
>Even if the data set was not generated properly, GRASS should somehow
>deal with that (as I cannot ask GLCF to redo the file).

GRASS coordinates in header files are corners, not center of cells, at
least that is how r.proj works. One half of the resolution is added to
get coordinates of the center of the corner cells. This is a good thing
(tm) because then we can avoid projection of nono points in some
projections, like 90N, but still define a location like 180W 90N 180E 90S.

I have some faint memory of DEM and GLCF being different in this
respect, DEM reported center of corner cells while GLCF used edges, but
I am not sure. Anyway, the user should check that the imported number of
columns/rows and resolution matches the original. Sometimes, if you do
the import with r.in.bin, half a resolution unit has to be
added/subtracted on each edge. r.in.gdal should get it right without
manual interference, but that probably depends on correct info in input
file header.

Thanks for your clarifying comments.

However, in this case it was r.in.gdal which caused problems:

http://grass.itc.it/pipermail/grass5/2005-April/017946.html
"r.in.gdal SRTM_GTOPO_u30_n090w020.tif out=srtm90_gtopo30
  ...
  D0/0: North: 90.0000000002
  WARNING: G_set_window(): Illegal latitude for North
"

Would be nice to accept such subtle deviations as shown above.

Using listgeo I get:

listgeo SRTM_GTOPO_u30_n090w020.tif
Geotiff_Information:
   Version: 1
   Key_Revision: 1.0
   Tagged_Information:
      ModelTiepointTag (2,3):
         0 0 0
         -20 90 0
      ModelPixelScaleTag (1,3):
         0.00833333377 0.00833333377 0
      End_Of_Tags.
   Keyed_Information:
      GTModelTypeGeoKey (Short,1): ModelTypeGeographic
      GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
      GeographicTypeGeoKey (Short,1): GCS_WGS_84
      GeogLinearUnitsGeoKey (Short,1): Linear_Meter
      GeogAngularUnitsGeoKey (Short,1): Linear_Meter
      End_Of_Keys.
   End_Of_Geotiff.

GCS: 4326/WGS 84
Datum: 6326/World Geodetic System 1984
Ellipsoid: 7030/WGS 84 (6378137.00,6356752.31)
Prime Meridian: 8901/Greenwich (0.000000/ 0d 0' 0.00"E)

Corner Coordinates:
Upper Left ( 20d 0' 0.00"W, 90d 0' 0.00"N)
Lower Left ( 20d 0' 0.00"W, 39d59'59.99"N)
Upper Right ( 20d 0' 0.01"E, 90d 0' 0.00"N)
Lower Right ( 20d 0' 0.01"E, 39d59'59.99"N)
Center ( 0d 0' 0.00"E, 64d60' 0.00"N)

As we see GDAL reports the contents of the file (which seems
to have some problem).

What about a flag in r.in.gdal to accept for epsilon deviation?
Then the user knows what he/she is doing.

Markus

14

On Mon, 11 Apr 2005 17:05:57 +0200
Markus Neteler <neteler@itc.it> wrote:

However, in this case it was r.in.gdal which caused problems:

http://grass.itc.it/pipermail/grass5/2005-April/017946.html
"r.in.gdal SRTM_GTOPO_u30_n090w020.tif out=srtm90_gtopo30
  ...
  D0/0: North: 90.0000000002
  WARNING: G_set_window(): Illegal latitude for North
"

Would be nice to accept such subtle deviations as shown above.

..

Corner Coordinates:
Upper Left ( 20d 0' 0.00"W, 90d 0' 0.00"N)

Note that is 22e-06 meters. We can make the test 0.1mm above 90deg lat,
0.0001 / (60*1852.) = 8.9993e-10

but will doing this allow weird segfaults with e.g. *.proj etc. down the
road? I think it might; better to fix the GDAL casting?

Hamish

15

> > I have noticed in the past that v.out.ogr to a shapefile does
> > similar strange off-by 1e-9 type errors to double precision columns
> > in the data.
>
> Mhh, a systematic problem?

No I think it is probably separate but wanted to note. A double cast to
float and then back to double introducing some noise?

example:

G6> echo "1.0|2.0|10.2" | v.in.ascii out=ogr_testG
...
column: 1 type: double
column: 2 type: double
column: 3 type: double
...

G6> dbview -bt $MAPSET/dbf/ogr_testG.dbf
1:1.000000:2.000000:10.200000:

[db.describe output below]

G6> v.out.ogr input=ogr_testG type=point dsn=. olayer=ogr_testS

G6> dbview -bt ogr_testS.dbf
1:1.000000000000000:2.000000000000000:10.199999999999999:

Hamish

G6> db.describe ogr_testG
table:ogr_testG
description:
insert:yes
delete:yes
ncols:4

column:cat
description:
type:INTEGER
len:11
scale:0
precision:10
default:0
nullok:yes
select:yes
update:yes

column:dbl_1
description:
type:DOUBLE PRECISION
len:20
scale:6
precision:18
default:0
nullok:yes
select:yes
update:yes

column:dbl_2
description:
type:DOUBLE PRECISION
len:20
scale:6
precision:18
default:0
nullok:yes
select:yes
update:yes

column:dbl_3
description:
type:DOUBLE PRECISION
len:20
scale:6
precision:18
default:0
nullok:yes
select:yes
update:yes

16

On Tue, 2005-04-12 at 18:37 +1200, Hamish wrote:

On Mon, 11 Apr 2005 17:05:57 +0200
Markus Neteler <neteler@itc.it> wrote:

>
> However, in this case it was r.in.gdal which caused problems:
>
> http://grass.itc.it/pipermail/grass5/2005-April/017946.html
> "r.in.gdal SRTM_GTOPO_u30_n090w020.tif out=srtm90_gtopo30
> ...
> D0/0: North: 90.0000000002
> WARNING: G_set_window(): Illegal latitude for North
> "
>
> Would be nice to accept such subtle deviations as shown above.
..
> Corner Coordinates:
> Upper Left ( 20d 0' 0.00"W, 90d 0' 0.00"N)

Note that is 22e-06 meters. We can make the test 0.1mm above 90deg lat,
0.0001 / (60*1852.) = 8.9993e-10

but will doing this allow weird segfaults with e.g. *.proj etc. down the
road? I think it might; better to fix the GDAL casting?

I did some comparisons between GRASS, listgeo and GDAL yesterday with
the above image.

In this case, GRASS is reporting wrong values (however small the
difference is). The projection handling needs to be reviewed.

--
Brad Douglas <rez@touchofmadness.com>