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

[Markus Neteler]

(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