[GRASSLIST:5704] Re: r.in.gdal problem?

[Glynn Clements]

Markus Neteler wrote:

> > The GRASS region refers to cell corners, not centres; it seems that
> > r.in.gdal is misinterpreting the region information in the grid file.
> > 
> > You can fix the resulting map manually using r.region, e.g.
> > 
> > 	r.region map=... n=77 s=-90 ...
> 
> Is this a GDAL bug or a GRASS bug? Does this mean that r.in.gdal
> is always shifting the data (were a nightmare!)?

I would think that it has to be a bug in either GDAL or r.in.gdal. If
the bug was in G_set_window(), I think that we would have noticed
before now ;)

r.in.gdal calls GDALGetGeoTransform(); if that succeeds (i.e. the data
is georeferenced), the returned values are used to initialise the
region. r.in.gdal doesn't add half-cell offsets.

If GDAL returns pixel centres, then r.in.gdal ought to be adding
half-cell offsets. However, gdal_datamodel.html says:

	The affine transform consists of six coefficients returned by
	GDALDataset::GetGeoTransform() which map pixel/line
	coordinates into georeferenced space using the following
	relationship:
	
	    Xgeo = GT(0) + Xpixel*GT(1) + Yline*GT(2)
	    Ygeo = GT(3) + Xpixel*GT(4) + Yline*GT(5)
	
	In case of north up images, the GT(2) and GT(4) coefficients
	are zero, and the GT(1) is pixel width, and GT(5) is pixel
	height. The (GT(0),GT(3)) position is the TOP LEFT CORNER of
	the top left pixel of the raster.

[Emphasis mine.]

So, it looks like a bug in GDAL.

AIGReadBounds() in frmts/aigrid/gridlib.c reads the values from the
dblbnd.adf file directly into psInfo->df{LL,UR}{X,Y}:

    VSIFRead( adfBound, 1, 32, fp );

	...
    
    psInfo->dfLLX = adfBound[0];
    psInfo->dfLLY = adfBound[1];
    psInfo->dfURX = adfBound[2];
    psInfo->dfURY = adfBound[3];

AIGDataset::GetGeoTransform() in frmts/aigrid/aigdataset.cpp adds the
half-cell offsets:

    padfTransform[0] = psInfo->dfLLX - psInfo->dfCellSizeX*0.5;
    padfTransform[1] = psInfo->dfCellSizeX;
    padfTransform[2] = 0;

    padfTransform[3] = psInfo->dfURY + psInfo->dfCellSizeY*0.5;
    padfTransform[4] = 0;
    padfTransform[5] = -psInfo->dfCellSizeY;

Presumably it shouldn't be doing this.

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