GeoTIFF: LANDSAT7 from GLCF projection info broken

[Markus Neteler]

Frank,

I am currently trying to read in
LANDSAT 7 from GLCF Maryland. This was successful some
time ago, maybe they have data in different formats.

From the site
ftp://ftp.glcf.umiacs.umd.edu/glcf/Landsat/WRS2/p195/r028/L71195028_02820010730.ETM-USGS/
e.g. (smallest):
L71195028_02820010730_B61.tif

reports:
gdalinfo L71195028_02820010730_B61.tif
Driver: GTiff/GeoTIFF
Size is 3410, 3157
Coordinate System is `'
GCP Projection = PROJCS["WGS 84 / UTM zone
32N",GEOGCS["WGS84",DATUM["WGS_1984",SPHEROID["WGS84",6378137,298.2572235629972,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",9],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32632"]]
GCP[  0]: Id=1, Info=
          (0,0) -> (297036,5.20918e+06,0)
Metadata:
  TIFFTAG_IMAGEDESCRIPTION={
  bandList =
  [
    6;
  ]
}
  TIFFTAG_DATETIME=2002:07:03 14:53:01

Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 3157.0)
Upper Right ( 3410.0,    0.0)
Lower Right ( 3410.0, 3157.0)
Center      ( 1705.0, 1578.5)
Band 1 Block=3410x307 Type=Byte, ColorInterp=Gray

While 'listgeo' tells me:
listgeo L71195028_02820010730_B61.tif
L71195028_02820010730_B61.tif: Warning, unknown field with tag 34595 (0x8723)
+encountered.
Geotiff_Information:
   Version: 1
   Key_Revision: 1.0
   Tagged_Information:
      ModelTiepointTag (2,3):
         0                0                0
         297036.055       5209182.38       0
      ModelTransformationTag (4,4):
         55.5416851       -12.8109804      0                297036.055
         -12.8109804      -55.5416851      0                5209182.38
         0                0                0                0
         0                0                0                1
      End_Of_Tags.
   Keyed_Information:
      GTModelTypeGeoKey (Short,1): ModelTypeProjected
      GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
      GTCitationGeoKey (Ascii,25): "Corrected Satellite Data"
      ProjectedCSTypeGeoKey (Short,1): PCS_WGS84_UTM_zone_32N
      PCSCitationGeoKey (Ascii,12): "UTM / WGS84"
      End_Of_Keys.
   End_Of_Geotiff.

PCS = 32632 (WGS 84 / UTM zone 32N)
Projection = 16032 (UTM zone 32N)
Projection Method: CT_TransverseMercator
   ProjNatOriginLatGeoKey: 0.000000 (  0d 0' 0.00"N)
   ProjNatOriginLongGeoKey: 9.000000 (  9d 0' 0.00"E)
   ProjScaleAtNatOriginGeoKey: 0.999600
   ProjFalseEastingGeoKey: 500000.000000 m
   ProjFalseNorthingGeoKey: 0.000000 m
GCS: 4326/WGS 84
Datum: 6326/World Geodetic System 1984
Ellipsoid: 7030/WGS 84 (6378137.00,6356752.31)
Prime Meridian: 8901/Greenwich (0.000000/  0d 0' 0.00"E)
Projection Linear Units: 9001/metre (1.000000m)

Corner Coordinates:
Upper Left    ( 297036.055,5209182.378)  (  6d19'48.01"E, 47d 0'18.08"N)
Lower Left    ( 256591.790,5033837.278)  (  5d53'20.64"E, 45d24'56.44"N)
Upper Right   ( 486433.201,5165496.935)  (  8d49'21.82"E, 46d38'34.43"N)
Lower Right   ( 445988.936,4990151.835)  (  8d18'50.27"E, 45d 3'45.96"N)
Center        ( 371512.496,5099667.107)  (  7d20'22.15"E, 46d 2'18.65"N)

The LANDSAT is covering Torino, North-West Italy. So the coordinates
printed by 'listgeo' sound reasonable.

Maybe you could look into this.

Kind regards

 Markus Neteler

[warmerdam]

Markus,

The basic problem with this dataset is that it provides a single tiepoint
*and* a transformation matrix.  GDAL's decision was to treat the tiepoint
as a GCP and to let that be the definitative georeferencing. 

I have modified the code to attach GCPs from tiepoints in addition to 
using the transformation matrix to define a geotransform with the additional
caveat that single tiepoints are ignored, and not transferred to the GCP
list. 

The changes are committed in CVS now. 

[neteler@…]

Frank,

thanks for the quick fix. Coordinates are reported now.
Unfortunately the resulting map appears to be distored:

http://mpa.itc.it/markus/gdf/ct/landsat7_globe30_vmap0.jpg
- Aspect from GLOBE30, reprojected to UTM of LANDSAT7
- Vector roads and rivers from VMAP0, reprojected to UTM of LANDSAT7
  (both align well)
- RGB composite of 1,2,3 superimposed.

I was using r.in.gdal (latest with today's GDAL-CVS) to import, it also
generated the GRASS location. Worked fine for other maps.

The large lake at North-West is shifted (warped?) by 20km.

If GDAL is right (simply applying the coordinates) it must be an
issue of the data vendor (GCLF University of Maryland/NASA/USGS).

Strange problem!

Best regards

 Markus Neteler

[warmerdam]

Markus,

The file will have a non-northup geotransform.  As far as I know, r.in.gdal
just ignores the rotational coefficients ... certainly there isn't a way of
representing rotated images in GRASS as far as I know.  

I would suggest you "launder" the files with gdalwarp before importing them
with r.in.gdal.  gdalwarp with default parameters should resample the files
into a north-up aspect.  

Ideally r.in.gdal should check for rotational coefficients, and at least
report the issue if they are non-zero. I am re-opening this bug report and
assigning it to Andrey who I have asked to take on r.in.gdal maintenance.
However, Andrey has a substantial contract of his own to work on now, and it
may be a while till he can look into this issue. 

[neteler@…]

Frank, Andrey,

the 'gdalwarp' did the trick. Now it looks very good.

I have already patched r.in.gdal in CVS:

+        if (adfGeoTransform[2] != 0.0 || adfGeoTransform[4] != 0.0)
+          G_fatal_error("Input map is rotated - cannot import. You may use
'gdalwarp' to transform the map to North-up.");
+

Hope that's fine with you.

 Markus

[warmerdam]

Markus, 

Excellent!  I will close this bug report now then. 

Best regards,