[GRASSLIST:5695] r.in.gdal problem?

David_Orme · February 28, 2003, 5:51pm

Hi,

I'm using the r.in.gdal module to import ARCINFO GRID format raster files and I seem to be getting a shift in the extent of the data on import:

The files contain 1 by 1 degree resolution raster data. In Arc, the NS extent of the raster data is:

bottom edge of southernmost grid square = 90S
topmost edge of northernmost square = 77N

The raster therefore contains 167 rows of data.

When I import this, the file extent is as follows:

GRASS 5.0.0 > r.info temp
+----------------------------------------------------------------------------+
| Layer: temp Date: Fri Feb 28 10:14:32 2003 |
| Mapset: dorme Login of Creator: dorme |
| Location: global |
| DataBase: /Users/dorme/grass |
| Title: ( temp ) |
|----------------------------------------------------------------------------|
| |
| Type of Map: raster Number of Categories: 1 |
| Data Type: CELL |
| Rows: 167 |
| Columns: 557 |
| Total Cells: 93019 |
| Projection: UTM (zone 13) |
| N: 77.5 S: -89.5 Res: 1 |
| E: 330.5 W: -226.5 Res: 1 |
| Range of data: min = 0 max = 1 |
| |
| Data Source: |
| |
| |
| |
| Data Description: |
| generated by r.in.gdal |
| |
| |
+----------------------------------------------------------------------------+

The number of rows (and columns) are correct but the extent is peculiar - I'm guessing that r.info returns the centre points of the grid cells but if this is the case then the centre of the northernmost cell should be at 76.5N not 77.5N. This also makes the reported southern extent odd since the 89.5 S implies that there is still a gridsquare that has 90S as a bottom edge. Using d.where to inspect the location visually shows that the top of the northernmost cell in grass is now at 78N and the bottom of the lowermost cell is at 89S.

Something doesn't add up - I'm expecting to find that it may be me but can anyone clear this up?

System details:
Grass 5.0.0 running under Mac OSX (10.2.4) within Apple's X11 (v.2)
Many thanks,
David

----------------------------------------
Dr. David Orme

Department of Biological Sciences,
Imperial College London,
Silwood Park Campus,
Ascot, Berkshire.
SL5 7PY

Phone: +44 (0)20 759 42358
e-mail: d.orme@imperial.ac.uk

Glynn_Clements · February 28, 2003, 8:39pm

David Orme wrote:

I'm using the r.in.gdal module to import ARCINFO GRID format raster
files and I seem to be getting a shift in the extent of the data on
import:

The files contain 1 by 1 degree resolution raster data. In Arc, the NS
extent of the raster data is:

bottom edge of southernmost grid square = 90S
topmost edge of northernmost square = 77N

The raster therefore contains 167 rows of data.

When I import this, the file extent is as follows:

> | Projection: UTM (zone 13) |
> | N: 77.5 S: -89.5 Res: 1 |
> | E: 330.5 W: -226.5 Res: 1 |

The number of rows (and columns) are correct but the extent is peculiar
- I'm guessing that r.info returns the centre points of the grid cells
but if this is the case then the centre of the northernmost cell should
be at 76.5N not 77.5N. This also makes the reported southern extent
odd since the 89.5 S implies that there is still a gridsquare that has
90S as a bottom edge. Using d.where to inspect the location visually
shows that the top of the northernmost cell in grass is now at 78N and
the bottom of the lowermost cell is at 89S.

Something doesn't add up - I'm expecting to find that it may be me but
can anyone clear this up?

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 ...

However, you appear to have imported the data into a UTM location,
when the data is actually lat/lon.

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

Markus_Neteler1 · March 1, 2003, 3:41pm

On Fri, Feb 28, 2003 at 08:39:35PM +0000, Glynn Clements wrote:

David Orme wrote:

> I'm using the r.in.gdal module to import ARCINFO GRID format raster
> files and I seem to be getting a shift in the extent of the data on
> import:
>
> The files contain 1 by 1 degree resolution raster data. In Arc, the NS
> extent of the raster data is:
>
> bottom edge of southernmost grid square = 90S
> topmost edge of northernmost square = 77N
>
> The raster therefore contains 167 rows of data.
>
> When I import this, the file extent is as follows:

> > | Projection: UTM (zone 13) |
> > | N: 77.5 S: -89.5 Res: 1 |
> > | E: 330.5 W: -226.5 Res: 1 |

> The number of rows (and columns) are correct but the extent is peculiar
> - I'm guessing that r.info returns the centre points of the grid cells
> but if this is the case then the centre of the northernmost cell should
> be at 76.5N not 77.5N. This also makes the reported southern extent
> odd since the 89.5 S implies that there is still a gridsquare that has
> 90S as a bottom edge. Using d.where to inspect the location visually
> shows that the top of the northernmost cell in grass is now at 78N and
> the bottom of the lowermost cell is at 89S.
>
> Something doesn't add up - I'm expecting to find that it may be me but
> can anyone clear this up?

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!)?

Markus

Glynn_Clements · March 1, 2003, 8:19pm

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@virgin.net>

John_Doucette · March 2, 2003, 10:38am

Hi everyone:

I am trying to set up a SummaSketch II Professional 12x18" digitizing tablet
to use with Grass 5.0.
I am running Redhat 71. Linux.
I have tried just plugging the tablet into the back of my computer and
choosing Cal Comp digitizer in the v.digit setup. I get an error message
suggesting I look in "/opt/..." for the existance of a "lock"
file/directory.

I don't have a clue how to proceed to configure the tablet. Can anyone
help?

John Doucette

Roy_Sanderson · March 3, 2003, 2:22pm

Dear John

I'm afraid that I don't have a simple answer as to how to configure your
digitizing tablet, as I have a Drawingboard III digitizer that isn't
talking to Grass5 in RH7.3. Presumably this is not straightforward in
Grass for some digitizer models or linux setups. However, if it is any
help to you, the key things you need to check are:

1) that your /etc/inittab file has a line referring to the correct tty, so
com port 1 is ttyS0, com 2 is ttyS1 e.g.
S1:235:off:/sbin/mingetty ttyS1 9600
would point to com 2. The man page on mingetty recommends using mgetty for
serial lines (although it made no difference to my setup).

2) that your /grass5/etc/digcap file is pointing to the correct digitizer
setup file as described in the /grass5/etc/dicap/digitizers directory. The
text descriptors for the available digitizers that appear when you run
v.digit are also in the digcap file.

3) that your digitizer is correctly configured to match one of the
digitizer setup files. From memory I seem to recall that the SummaSketch
uses a series of dip-switches on the underside of the digitizer to control
its baud rate, data mode etc.

Good luck.
Best wishes
Roy

At 02:38 02/03/03 -0800, you wrote:

Hi everyone:

I am trying to set up a SummaSketch II Professional 12x18" digitizing tablet
to use with Grass 5.0.
I am running Redhat 71. Linux.
I have tried just plugging the tablet into the back of my computer and
choosing Cal Comp digitizer in the v.digit setup. I get an error message
suggesting I look in "/opt/..." for the existance of a "lock"
file/directory.

I don't have a clue how to proceed to configure the tablet. Can anyone
help?

John Doucette

----------------------------------------------------------------------------
Roy Sanderson
Centre for Life Sciences Modelling
Porter Building
University of Newcastle
Newcastle upon Tyne
NE1 7RU
United Kingdom

Tel: +44 191 222 7789

r.a.sanderson@newcastle.ac.uk
http://www.ncl.ac.uk/clsm

----------------------------------------------------------------------------