[GRASS5] r.in.gdal and negative values

GRASS GIS GRASS GIS developer
,
1

re. http://grass.itc.it/pipermail/grass5/2005-January/017029.html

Currently r.in.gdal imports unprojected data in the negative, i.e. 0,0
is the top right of the image.

(unlike r.in.tiff|png in GRASS 5.4 which makes 0,0 bottom left)

This can be "fixed" with r.region by the user after import.

some comments from the source code:

/* TODO: most unreferenced formats are read in with negative
* coordinates - desired??
*/

[...]

        /* use negative XY coordinates per default for unprojected data
           but not for all formats... (MN: I don't like it at all) */
        /* for hDriver names see gdal/frmts/gdalallregister.cpp */

        if ( l1bdriver ||
             ( strcmp(GDALGetDriverShortName(hDriver),"GTiff") ||
               strcmp(GDALGetDriverShortName(hDriver),"JPEG") ) == 0 )
        {
          /* e.g. L1B - NOAA/AVHRR data must be treated differently */
          /* define positive xy coordinate system to avoid GCPs confusion */
          cellhd.north = cellhd.rows;
          cellhd.south = 0.0;

[...]
        else
        {
          /* for all other unprojected data ... */
          /* define negative xy coordinate system to avoid GCPs confusion */
          cellhd.north = 0.0;
          cellhd.south = (-1) * cellhd.rows;

I don't understand the logic or intention here: who's supposed to be
getting set to negative and why?

FWIW the above strcmp() part of the test will always fail.

i.e. strcmp() returns 0 if the strings match. It will never match both a
GeoTIFF and a JPEG, so one side or the other of the "OR" will be always
be non-zero and the ==0 part will always be false. So it really reduces
to if(l1bdriver) {make positive} else {make negative}.

(l1bdriver is AVHRR; http://www.gdal.org/frmt_l1b.html
  Will this data ever be without GCPs & so touch this part of the code?!)

[And PNG should be doing the same thing as TIFF & JPEG.]

So what is the intent of the test?

Are GeoTIFFs and JPEGs supposed to be ending up positive or negative?
  (positive please)

Who's getting what confused with GCPs? Does that situation still exist?

I prefer to make everything positive as it's less confusing.

gdalinfo reports this on an unprojected TIFF:
Driver: GTiff/GeoTIFF
Size is 1309, 4920
...
Corner Coordinates:
Upper Left ( 0.0, 0.0)
Lower Left ( 0.0, 4920.0)
Upper Right ( 1309.0, 0.0)
Lower Right ( 1309.0, 4920.0)
Center ( 654.5, 2460.0)

So in GDAL 0,0 is top left (image coordinates (or satellite scanlines)
as opposed to cartesian coordinates)

GRASS XY uses cartesian coordinates so 0,0 should be bottom left if we
want to keep map units positive. (external support for idea: this is the
reason for having false eastings and northings)

So is this the "GCP confusion" the code speaks of? Solution is to make
it neither? :confused: As we are GRASS, matching GRASS XY conventions seems
like logical thing to do.

I would like to use i.point's $GROUP/POINTS file as a template for
setting GCPs in a GeoTIFF directly with gdal_translate or geotifcp,
but it seems like I'll have to fuddle the i.points results around
slightly to do so? (flip y-axis) If everything is positive it is
easy enough to write a script to form a "gdal_translate -gcp"
command string from i.points output.

It seems it is trivial to make everything positive in the r.in.gdal
code; can anyone justify a reason not to?

Aka any reason I should not just go ahead and do this in 6.1-cvs?

thanks,
Hamish

2

Hamish schrieb:

It seems it is trivial to make everything positive in the r.in.gdal code; can anyone justify a reason not to?

Aka any reason I should not just go ahead and do this in 6.1-cvs?

This is one of the excrescent obstacles when grass-novices try to go their fist steps!

In particular because you will neither get the info that your import was successful with negative coordinates nor you see anything when calling display function.

So it would be very helpful if you could fix it.

Mit freundlichen Grüßen / With kindest regards

Stefan Paulick

http://www.urbeli.com
mailto://stefan.paulick@urbeli.com
/*----------------------*/

3

On Tue, May 31, 2005 at 05:06:26PM +1200, Hamish wrote:

re. http://grass.itc.it/pipermail/grass5/2005-January/017029.html

Currently r.in.gdal imports unprojected data in the negative, i.e. 0,0
is the top right of the image.

(unlike r.in.tiff|png in GRASS 5.4 which makes 0,0 bottom left)

This can be "fixed" with r.region by the user after import.

some comments from the source code:

/* TODO: most unreferenced formats are read in with negative
* coordinates - desired??
*/

[...]

        /* use negative XY coordinates per default for unprojected data
           but not for all formats... (MN: I don't like it at all) */

... so you can see that we are two people not being in favor of
negative xy coordinates...

        /* for hDriver names see gdal/frmts/gdalallregister.cpp */

        if ( l1bdriver ||
             ( strcmp(GDALGetDriverShortName(hDriver),"GTiff") ||
               strcmp(GDALGetDriverShortName(hDriver),"JPEG") ) == 0 )
        {
          /* e.g. L1B - NOAA/AVHRR data must be treated differently */
          /* define positive xy coordinate system to avoid GCPs confusion */
          cellhd.north = cellhd.rows;
          cellhd.south = 0.0;

[...]
        else
        {
          /* for all other unprojected data ... */
          /* define negative xy coordinate system to avoid GCPs confusion */
          cellhd.north = 0.0;
          cellhd.south = (-1) * cellhd.rows;

I don't understand the logic or intention here: who's supposed to be
getting set to negative and why?

As written: L1B data have a strange GCP organization. We should
probably change the test to refuse L1B data completely as
i.rectify is not suitable to recitfy AVHRR data.

The new gdalwarp/thin plate splines warper performs better (so
a new error message could suggest exactly this).

FWIW the above strcmp() part of the test will always fail.

Mhh, it worked at least on a couple of platforms.
Since I wrote it it could be certainly programmed in a better way.

i.e. strcmp() returns 0 if the strings match. It will never match both a
GeoTIFF and a JPEG, so one side or the other of the "OR" will be always
be non-zero and the ==0 part will always be false. So it really reduces
to if(l1bdriver) {make positive} else {make negative}.

This was probably intended (Frank wanted always negative coordinates which
didn't work for AVHRR). So I implemented the exception for L1B.

(l1bdriver is AVHRR; http://www.gdal.org/frmt_l1b.html
  Will this data ever be without GCPs & so touch this part of the code?!)

It's always with GCPs which are written into a POINTS file during
import.

[And PNG should be doing the same thing as TIFF & JPEG.]

Yes.

So what is the intent of the test?

Are GeoTIFFs and JPEGs supposed to be ending up positive or negative?
  (positive please)

Frank wanted negative (I don't know why - maybe there is no reason).
I suggest to change to positive values for all and to refuse L1B or,
at least, issue a warning.

Who's getting what confused with GCPs? Does that situation still exist?

I prefer to make everything positive as it's less confusing.

Agreed.

gdalinfo reports this on an unprojected TIFF:
Driver: GTiff/GeoTIFF
Size is 1309, 4920
...
Corner Coordinates:
Upper Left ( 0.0, 0.0)
Lower Left ( 0.0, 4920.0)
Upper Right ( 1309.0, 0.0)
Lower Right ( 1309.0, 4920.0)
Center ( 654.5, 2460.0)

So in GDAL 0,0 is top left (image coordinates (or satellite scanlines)
as opposed to cartesian coordinates)

Note that GDAL uses *corners* of pixels for coordinates.

GRASS XY uses cartesian coordinates so 0,0 should be bottom left if we
want to keep map units positive. (external support for idea: this is the
reason for having false eastings and northings)

Well, false eastings and northings do not appear in XY locations AFAIK.

So is this the "GCP confusion" the code speaks of? Solution is to make
it neither? :confused: As we are GRASS, matching GRASS XY conventions seems
like logical thing to do.

The story is:
- Frank wrote the module with negative values in the XY case
- I added special treatment for AVHRR to make it work. Then I
  discovered that polynomials perform bad for AVHRR. Finally
  I organized partial fundings for the thin spline warper in GDAL
  which works way better.

I would like to use i.point's $GROUP/POINTS file as a template for
setting GCPs in a GeoTIFF directly with gdal_translate or geotifcp,
but it seems like I'll have to fuddle the i.points results around
slightly to do so? (flip y-axis) If everything is positive it is
easy enough to write a script to form a "gdal_translate -gcp"
command string from i.points output.

It seems it is trivial to make everything positive in the r.in.gdal
code; can anyone justify a reason not to?

Aka any reason I should not just go ahead and do this in 6.1-cvs?

Please go ahead...!

Markus

4

> re. http://grass.itc.it/pipermail/grass5/2005-January/017029.html
>
>
> Currently r.in.gdal imports unprojected data in the negative, i.e.
> 0,0 is the top right of the image.
>
> (unlike r.in.tiff|png in GRASS 5.4 which makes 0,0 bottom left)
>
> This can be "fixed" with r.region by the user after import.

...

> I don't understand the logic or intention here: who's supposed to be
>
> getting set to negative and why?

As written: L1B data have a strange GCP organization. We should
probably change the test to refuse L1B data completely as
i.rectify is not suitable to recitfy AVHRR data.

The new gdalwarp/thin plate splines warper performs better (so
a new error message could suggest exactly this).

...

Frank wanted negative (I don't know why - maybe there is no reason).
I suggest to change to positive values for all and to refuse L1B or,
at least, issue a warning.

I think a warning is best. GRASS XY may still be used as a data-
exploration or translation software. In future i.rectify may be updated
with more methods. As a general rule, I don't like to enforce that
something is not possible just because I don't see why it is important
or how to do it myself.

The story is:
- Frank wrote the module with negative values in the XY case

I do remember a thread about this on the mailing list some time ago but
can't seem to find it in the archives. Bad choice of search words I
guess. Any hints? If we can't find Frank's reasoning, perhaps we should
wait a few days until he is back on-line. And then do it anyway :slight_smile:

Found this:
https://intevation.de/rt/webrt?serial_num=805

- I added special treatment for AVHRR to make it work. Then I
  discovered that polynomials perform bad for AVHRR. Finally
  I organized partial fundings for the thin spline warper in GDAL
  which works way better.

I plan on testing out the GDAL splines soon, will let you know how it
goes.

cheers,
Hamish

5

I suggest to change to positive values for all and to refuse L1B or,
at least, issue a warning.

Done.

L1B allowed, but with warning to use 'gdalwarp -tps' instead of i.rectify.

Hamish