Coordinate system not reported for ERDAS .IMG file

[Markus Neteler]

Hi Frank,

at time we face a bit troubles with .IMG files.
It seems that 'gdalinfo' (latest CVS) does not show the 
correct coordinate system.

Apparently there is a problem in GDAL because ArcGIS
reads the .img files perfectly. I'm sending you two
.gif  to show what I mean. A colleague opened the two
georeferenced images in ArcGIS and put over it an
Europe map. Both images are ch2 AVHRR. The green one
used a 1st order transformation and the red used a 2nd
order.Notice the close-up image on the Spain/France border,
how the ocean contour fits nicely on the 2nd order
transformation. So it appears that ERDAS is working
fine. These ARCGIS snapshot images are stored here:

gdal.velocet.ca/pub/incoming/
 erdas_arcgis_avhrrimages.gif
 erdas_arcgis_avhrrspain_close.gif

Trying the same with GRASS/GDAL it fails because only
xy coordinates are reported.
The original erdas files (generated with Imagine) are
gdal.velocet.ca/pub/incoming/
 erdasavhrr_2nd.img
 erdasavhrr_2nd.rrd

Here is what I get:
gdalinfo erdasavhrr_2nd.img
Driver: HFA/Erdas Imagine Images (.img)
Size is 2708, 261
Coordinate System is:
GEOGCS["WGS 72",
    DATUM["WGS 72",
        SPHEROID["WGS 72",6378135,298.259998589995],
        TOWGS84[0,0,4.5,0,0,-2.686e-06,2.2e-07]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]]
Origin = (0.000000,0.000000)
Pixel Size = (1.000000,1.000000)
Corner Coordinates:
Upper Left  (   0.0000000,   0.0000000) (  0d 0'0.00"E,  0d 0'0.00"N)
Lower Left  (       0.000,     261.000) (  0d 0'0.00"E,261d 0'0.00"N)
Upper Right (    2708.000,       0.000) (2708d 0'0.00"E,  0d 0'0.00"N)
Lower Right (    2708.000,     261.000) (2708d 0'0.00"E,261d 0'0.00"N)
Center      (    1354.000,     130.500) (1354d 0'0.00"E,130d30'0.00"N)
Band 1 Block=64x64 Type=Float32, ColorInterp=Undefined
  Overviews: 675x63
Band 2 Block=64x64 Type=Float32, ColorInterp=Undefined
  Overviews: 675x63
Band 3 Block=64x64 Type=Float32, ColorInterp=Undefined
  Overviews: 675x63
Band 4 Block=64x64 Type=Float32, ColorInterp=Undefined
  Overviews: 675x63
Band 5 Block=64x64 Type=Float32, ColorInterp=Undefined
  Overviews: 675x63

It seems that a mixture of xy and Lat/Long coodinates is
reported.

Perhaps it is an user error, but since ArcGIS accepts the
erdasavhrr_2nd.img file, GDAL may also accept it...

Thanks in advance,

 Markus

[warmerdam]

Markus,

I have looked quickly at the file, and the problem is that the imagine
file holds the pixel to georeferenced transformation as special polynomial
transform structures.  This does not map cleanly into the GDAL concept of 
georeferencing mechanisms, though I suppose I could use the polynomials to
generate a grid of GCPs.  I will take this as a low priority enhancement 
request. 

I have also obserserved that "hfatest", my imagine dumper, fails on the
file.  It crashes dumping the polynomial information apparently because 
the code for extracting values doesn't properly understand how to skip
"BASEDATA".  This is left as a higher priority bug to fix, though it
doesn't affect you directly.

    MapToPixelXForm(Exfr_GenericXFormHeader) @ 1930 + 27 @ 2058
    + titleList = 
    +     string = `Polynomial'

      XForm0(GM_PolyPair) @ 2085 + 432 @ 17633376
      + forward = 
      +     order = 2
      +     numdimtransform = 2
      +     numdimpolynomial = 2
      +     termcount = 6
      +     exponentlist[0] = 0
      +     exponentlist[1] = 0
      +     exponentlist[2] = 1
      +     exponentlist[3] = 0
      +     exponentlist[4] = 0
      +     exponentlist[5] = 1
      +     exponentlist[6] = 2
      +     exponentlist[7] = 0
      +     exponentlist[8] = 1
      +     exponentlist[9] = 1
      +     exponentlist[10] = 0
      +     exponentlist[11] = 2
      +     polycoefmtx = (basedata)
      +     polycoefvector[0] = (basedata)
      +     polycoefvector[1] = (basedata)
      + reverse[0] = 
      +     order = -1067575556
      +     numdimtransform = 1745729341
      +     numdimpolynomial = -1069140783
      +     termcount = 242191230
Segmentation fault

HFAType Efga_Polynomial/-1 bytes
    ULONG                 order[1];
    ULONG                 numdimtransform[1];
    ULONG                 numdimpolynomial[1];
    ULONG                 termcount[1];
    ULONG               p exponentlist[0];
    BASEDATA            * polycoefmtx[1];
    BASEDATA            * polycoefvector[1];

[neteler@…]

Frank,

new LANDSAT TM5 IMG data are arrived... gdalinfo can partially get the projection:

gdalinfo 174073mssres11July2001.img
Driver: HFA/Erdas Imagine Images (.img)
Size is 8399, 7360
Coordinate System is:
PROJCS["UTM",
    GEOGCS["Undefined",
        DATUM["Undefined",
            SPHEROID["WGS 84",6378137,298.2572235629972]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    UNIT["meters",1],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-183],         <--- should be different
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",10000000]]
Origin = (692047.000000,8030859.000000)
Pixel Size = (30.000000,-30.000000)
Corner Coordinates:
Upper Left  (  692047.000, 8030859.000) (178d48'42.66"E, 17d48'5.01"S)
Lower Left  (  692047.000, 7810059.000) (178d49'59.98"E, 19d47'45.38"S)
Upper Right (  944017.000, 8030859.000) (178d48'51.45"W, 17d45'53.99"S)
Lower Right (  944017.000, 7810059.000) (178d45'53.35"W, 19d45'18.56"S)
Center      (  818032.000, 7920459.000) (179d58'59.93"W, 18d46'59.65"S)
Band 1 Block=64x64 Type=Byte, ColorInterp=Undefined
Band 2 Block=64x64 Type=Byte, ColorInterp=Undefined
Band 3 Block=64x64 Type=Byte, ColorInterp=Undefined
Band 4 Block=64x64 Type=Byte, ColorInterp=Undefined
Band 5 Block=64x64 Type=Byte, ColorInterp=Undefined
Band 6 Block=64x64 Type=Byte, ColorInterp=Undefined

The data are either UTM SOUTH zone 34 or 35 (Nord Botswana).
So the data should be located somewhere around

Center Longitude: 23:47:17.841853E [23.78829]
Center latitude:  18:32:35.276525S [-18.54313]

Maybe the 'gdalwarp' related improvements allow to get this
working (in future...).

The file is 370MB, so I don't upload it.

Best regards

 Markus

[warmerdam]

Markus, 

If you cd into gdal/frmts/hfa and do a "make hfadump" you should get a utility
which can report stuff about a .img file.  Then use "./hfadump -dt <yourfile>.img"
to get a dump of the internal tree structure of the file.  Append that output
to this bug report. 

Thanks,

[neteler@…]

Frank,

I have compiled hdatest and attached the output.
Here I found some tools from ERDAS, but couldn't compile on
Linux (there are SUN etc binaries as well, but I have only Linux).
ftp://ftp.erdas.com/pub/imgautil/

Best regards

 Markus

[neteler@…]

hfatest output of LANDSAT TM 5 data set

[neteler@…]

gdalinfo output of same IMG file

[warmerdam]

Markus, 

The "Projection" node is the issue here.  It looks like this:

Projection(Eprj_ProParameters) @ 359845386 + 221 @ 359845514
    + proType = EPRJ_INTERNAL
    + proNumber = 1
    + proExeName = (no values)
    + proName = `UTM'
    + proZone = 0
    + proParams[0] = 0.000000
    + proParams[1] = 0.000000
    + proParams[2] = 0.000000
    + proParams[3] = -1.000000
    + proParams[4] = 0.000000
    + proParams[5] = 0.000000

The problem is that the zone value isn't set properly. Does this UTM definition
work properly in Imagine or the erdas free viewer?  I don't see how it
possibly good be extracted properly.

Currently it is taking the zone number to be 0 which doesn't exist.  Note
that the "proParams[3] == -1" marks it as being in the southern hemisphere,
a fact which seems to be picked up properly.

I don't really see any way to fix this in GDAL.

Best regards,

[neteler@…]

Frank,

thanks for revisiting this problem.
I have tried on MS-Windows with ERDAS ViewFind, it looks
like "bad luck" for me. ViewFind reports the same as GDAL.

In fact the zone value isn't set properly.

I asked for the data provider - the data were prepared by an
official satellite data company.

So I should be able to fix that with gdal_translate and
EPSG 32735 (WGS 84 / UTM zone 35S).

gdal_translate -of GTiff -a_srs '+init=epsg:32735' -b 3
173072mssres18June2001.img 173072mssres18June2001.tif
....
Done.
And verified with other data - looks ok.

Thanks for providing gdal_translate - it seems that I am able
to fix these corrupted data now.
[you may close the bug]

Best regards

 Markus