[GRASS-user] r.in.gdal global raster incorrect extents

[Hamish]

> > Jamie Adams wrote:
> > > In my case, the source GeoTiff had a geotransform of:
> > > (-180.00000000012497, 0.0083333333333333297, 0.0, 90.0,
> > > 0.0, -0.0083333333333333297)
...
> My global 0.5 min grid is working, but I've discovered problems with
> other data sets

Hamish:
> > what is the problem besides a slightly-broken input file?

Jamie:
> I'm not sure why this is considered a broken input file.  Yes it has a
> slight shift to the west, but all of the cell centers lie between 180w
> 90n 180e 90s.  Yes, this amount of shift is trivial and can be easily
> fixed, but it could be considerably worse. Take for example some other
> raster data I have that is registered up to 180w, but has the pixel
> centers aligned to the grid (i think srtm is this way).  I do a
> r.in.gdal, g.region rast=in.tif, r.out.gdal and end up with a
> completely different registration.

 
> gdalinfo W180N10.tif
....
> Origin = (-180.000416666999996,10.000416666700000)
> Pixel Size = (0.000833333333333,-0.000833333333333)

note in this case the amount it exceeds 180 by is 1/2 the cell resolution.
i.e. the cells are centered on the grid lines. read about that here:
  http://grass.osgeo.org/wiki/GRASS_raster_semantics
  http://www.ngdc.noaa.gov/mgg/global/

Adjustments therefore need to be made so that coordinates refer to the
cell boundaries, not centers. For an idea how that is done, read through
the r.in.srtm script:
http://trac.osgeo.org/grass/browser/grass/trunk/scripts/r.in.srtm/r.in.srtm

So in this case, the SRTM tile is not broken, it is what it is and we
just have to process it correctly to compensate.


in your prior example though (-180.00000000012497, ...) the >180 error
is not due to grid vs. cell conventions (off by 1/2 a cell). It is much
more likely due to a loss of floating point precision in the program that
created it. Thus I call that one "broken". I don't think GDAL is to blame,
as for me GDAL reports the origin longitude cleanly. (see below).

What value gets stored in the $MAPSET/cellhd/ file for this image??



> r.in.gdal input=W180N10.tif output=W180N10
> 
> r.info -g W180N10
> 
> north=10:00:01.5N
> south=10:00:01.5S
> east=159:59:58.500001W
> west=179:59:58.499999E
> 
> g.region rast=W180N10
> 
> r.out.gdal input=W180N10 output=W180N10_v2.tif
> 
> gdalinfo W180N10_v2.tif
...
> Origin = (179.999583333055540,10.000416666666666)
> Pixel Size = (0.000833333333333,-0.000833333333333)

in the above the longitude is ok for at first, but then loses it for
the later digits: 179.999583333 055540.

To clean things up you might try the g.region -a flag. You may need to
double the resolution if that "rounds to whole numbers" and you want
the region to be "on the halves". [e.g. w=25 e=75 res=50 -a would "round"
outwards to w=0 e=100 res=50; it fails to notice the grid is already at
50m if the whole thing is offset. setting res=25 with -a avoids the 1/2
cell shift. you can set res back to 50 after the -a align step]


for a SRTMv2 tile I have here, it looks ok:
$ gdalinfo S46E170.hgt
....
Origin = (169.999583333333334,-44.999583333333334)
Pixel Size = (0.000833333333333,-0.000833333333333)

ie the origin longitude doesn't lose precision.  (GDAL 1.5.0)


make sense?

Hamish