Discrepancy between gdalinfo and QGIS raster metadata

[alobo]

At reading file test2.tif (use pseudocolor for visualization), which is available here ​http://sites.google.com/site/eospansite/dummy/test2.tif (sometimes the direct link does not work, in such a case point your browser to ​http://sites.google.com/site/eospansite/dummy and find the file, the 3rd one starting from the bottom)

the file info issued by gdalinfo (using QGIS plugin Raster/Information) is different than the one offered by QGIS under Raster Layer Properties/Metadata. The file was created by QGIS plugin Interpolation out of a vector points file in ED50 UTM31N projection. The fact is that the orientation in QGIS is correct, while R follows the info given by gdal and the image is reversed. Note that QGIS states the origin in the NW corner and then sets negative Y resolution. But this is not what gdalinfo reads.

Details:

gdalinfo test2.tif Driver: GTiff/GeoTIFF Files: test2.tif Size is 256, 256 Coordinate System is: PROJCS["ED50 / UTM zone 31N",

GEOGCS["ED50",

DATUM["European_Datum_1950",

SPHEROID["International 1924",6378388,297.000000000005,

AUTHORITY["EPSG","7022"]],

AUTHORITY["EPSG","6230"]],

PRIMEM["Greenwich",0], UNIT["degree",0.0174532925199433], AUTHORITY["EPSG","4230"]],

PROJECTIONTransverse_Mercator, PARAMETER["latitude_of_origin",0], PARAMETER["central_meridian",3], PARAMETER["scale_factor",0.9996], PARAMETER["false_easting",500000], PARAMETER["false_northing",0], UNIT["metre",1,

AUTHORITY["EPSG","9001"]],

AUTHORITY["EPSG","23031"]]

Origin = (424389.000000000000000,4635822.000000000000000) Pixel Size = (345.007812500000000,194.019531250000000) Metadata:

AREA_OR_POINT=Area

Image Structure Metadata:

INTERLEAVE=BAND

Corner Coordinates: Upper Left ( 424389.000, 4635822.000) ( 2d 5'20.10"E, 41d52'11.88"N) Lower Left ( 424389.000, 4685491.000) ( 2d 4'56.97"E, 42d19'2.06"N) Upper Right ( 512711.000, 4635822.000) ( 3d 9'11.42"E, 41d52'24.52"N) Lower Right ( 512711.000, 4685491.000) ( 3d 9'15.31"E, 42d19'14.90"N) Center ( 468550.000, 4660656.500) ( 2d37'10.89"E, 42d 5'47.83"N) Band 1 Block=256x4 Type=Float64, ColorInterp=Gray

while QGIS Raster Layer Properties/Metadata states: Layer Spatial Reference System: +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs Origin: 424389,4.68549e+06 Pixel Size: 279.887,-279.887

[mmassing]

Looking at the sample file you provided, I think the issue here is that the generated .tif is not north-up (i.e. the geotransform contains a non-zero rotational component or, in this case, a positive scale factor for the y-axis).

Such coordinate systems are perfectly valid, but some software may not be prepared to handle this general case.

QGis handles this type of files by reprojecting them on the fly (using the gdal VRT mechanism). Other applications which do not explicitly support a rotated/flipped coordinate system may require you to warp the input file into a north-up orientation (e.g. by using "gdalwarp test2.tif test2-nu.tif", or by wrapping it in an appropriate VRT file).

Manuel

[mmassing]

To summarize this more concisely, I think this is a bug in R (which fails to handle general coordinate systems), and probably can be worked around by reprojecting the problematic file to a north-up orientation.

[alobo]

I'll let Roger Bivand know about this thread, but beyond R, why is the info displayed by gdalinfo (not R) different than what is shown by the raster metadata in QGIS? Also, note that test2.tif was generated by the interpolation tool of QGIS from a points vector layer in UTM 31N ED50, why is this raster file "not north-up (i.e. the geotransform contains a non-zero rotational component or, in this case, a positive scale factor for the y-axis)." ?

Agus

[mmassing]

Displaying the VRT coordinate frame instead of the original raster coordinate frame can arguably be seen as a user interface issue, but I can't really comment on that.

Maybe Marco Hugentobler has any insight into why the orientation of the interpolation output has a positive y-axis scale (which, as already mentioned, is not strictly a bug, but might trip up some unprepared software - so if the y-axis scale can be easily adjusted by the interpolation plugin, this would probably be a good solution to this issue).

Manuel

[rbivand]

Replying to mmassing:

To summarize this more concisely, I think this is a bug in R (which fails to handle general coordinate systems), and probably can be worked around by reprojecting the problematic file to a north-up orientation.

No, in no way. The R rgdal package interfaces GDAL, and the GDAL driver sees the file in exactly the same way. If the origin is NW, but the resolution increment is positive, as here, any application using the information GDAL has will behave in the same way. No jumping through hoops here, the problem must be in QGIS and/or the plugin writing the file with inappropriate parameters. I'm not planning to make any changes in the R/GDAL interface to suit a very odd case. If the data were correctly saved and described, GDAL would see them correctly. Maybe the GDAL driver might accommodate oddities, but at out end, we can't make adjustments for one application.

[rbivand]

Replying to rbivand:

Replying to mmassing:

To summarize this more concisely, I think this is a bug in R (which fails to handle general coordinate systems), and probably can be worked around by reprojecting the problematic file to a north-up orientation.

No, in no way. The R rgdal package interfaces GDAL, and the GDAL driver sees the file in exactly the same way. If the origin is NW, but the resolution increment is positive, as here, any application using the information GDAL has will behave in the same way. No jumping through hoops here, the problem must be in QGIS and/or the plugin writing the file with inappropriate parameters. I'm not planning to make any changes in the R/GDAL interface to suit a very odd case. If the data were correctly saved and described, GDAL would see them correctly. Maybe the GDAL driver might accommodate oddities, but at out end, we can't make adjustments for one application.

In addition, it appears that QGIS is doing things in src/core/raster/qgsrasterlayer.cpp that perhaps use geoTransform values in interpretative ways, like switching Ymin and Ymax from those declared in the file. It would be better to declare them right initially, wouldn't it?

I cannot locate the offending interpolation tool or plugin to examine its source code. Could someone tell me where it is?

My understanding is that QGIS knows that the raster is saved the wrong way up, but uses VRT to display it correctly. In my opinion it simply should be stored the right way up, which would remove the need for jumping through hoops. I can consider adding a warning to rgdal raster reading to indicate that flip() should be used in the imported raster object, but there are limits to what is sensible.

[mhugent]

With which tool has this tiff (test2.tif) been generated? It cannot be with the interpolation C++ plugin, because this only writes ascii grid files as output (it is one of my intended mid-term improvements to save it in any supported GDAL format, such as GeoTiff).

Was it the GDAL Toos (python plugin)? Or the raster analysis (C++) plugin?

[rbivand]

Replying to mhugent:

With which tool has this tiff (test2.tif) been generated? It cannot be with the interpolation C++ plugin, because this only writes ascii grid files as output (it is one of my intended mid-term improvements to save it in any supported GDAL format, such as GeoTiff).

Was it the GDAL Toos (python plugin)? Or the raster analysis (C++) plugin?

This agrees with my analysis using 1.3.0 on F10. Agus said that the file came from 1.4.0 on Windows. If the interpolator tool/plugin only produces ARC ASCII files, then we do need to find out where the GTiff came from. The interpolator is not the culprit here.

The interpolator-generated Arc ASCII raster that I made and read into R using GDAL was imported correctly. I would, however, give a default square cell resolution in the interpolator menu, which would mean that the GDAL-specific DX/DY tags are not used (using a check box to revert to current behaviour). Further, it would be great to be able to set the power for IDW from the menu (maybe it is there, but I didn't see it).

[alobo]

Let me reproduce the process again this afternoon using the same machine. The interpolation was made with the interpolator plugin and I think that this raster was later written as geotif using translate in the gdal plugin, but let me make sure. I was wrong at guesing that the win version of the interpolator could have generated the geotif file.

Agus

[alobo]

ok, I think I've been able to clarify this:

1. Interpolation plugin generates an Arc Ascii raster, even if it's named test2.tif
2. This file has the correct orientation once imported in R through either readGDAL()

or raster()

3. But lacks projection information: no projection information is displayed by gdalinfo,

both through GDALinfo() and through Raster/Information in QGIS gdal plugin.

4. Then I used QGIS Raster/Assign projection to assign UTM31N ED50. The projection

is assigned and the file converted to geotiff (without asking you)

5. The new test2.tif shows correct projection information but is reversed once imported in R.

So the problem is in Raster/Assign projection. I still have to check if this happens with any Raster/ utility in QGIS, hope not.

Agus

[brushtyler]

I think there is not only one issue.

First issue:
Replying to alobo:

1. Interpolation plugin generates an Arc Ascii raster, even if it's named test2.tif

why Interpolation plugin create an Arc Ascii named .tif? Do you have opened a ticket about this?

Second (possible) issue:

4. Then I used QGIS Raster/Assign projection to assign UTM31N ED50. The projection

is assigned and the file converted to geotiff (without asking you)

the gdal Assign Projection gui says:

The output is: 
- new GeoTiff if input file is not GeoTiff
- overwritten if input is GeoTiff

It checks if is a GeoTiff by looking at the file extension, so it overwrite the input file.

Replying to alobo:

So the problem is in Raster/Assign projection. I still have to check if this happens with any Raster/ utility in QGIS, hope not.

Assign Projection use gdalwarp. Have you tried to change the projection with another gdal executable (from shell)?

[alobo]

I agree there is more than one issue, I was actually pointing 4 issues (2 and 3 are just two steps of the same problem). Probably the 4th one is a minor one, QGIS could just warn that a tif file will be written. Nevertheless, QGIS should follow its own consistent rules at writing files, and the best practice should always be asking the user for the output format. Regarding the 1st issue, I have to test again and, eventually, open a different ticket. This issue created a lot of confusion initially. 2&3 and 5 are still open problems, although I'd have to revisit the whole thing to make sure. Agus