Bug report #2444

Discrepancy between gdalinfo and QGIS raster metadata

Added by alobo - over 10 years ago. Updated over 1 year ago.

Status:Closed
Priority:Low
Assignee:-
Category:Rasters
Affected QGIS version:master Regression?:No
Operating System: Easy fix?:No
Pull Request or Patch supplied:No Resolution:end of life
Crashes QGIS or corrupts data:No Copied to github as #:12504

Description

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,
AUTHORITYEPSG""7022,
AUTHORITYEPSG""6230,
PRIMEMGreenwich,
UNITdegree,
AUTHORITYEPSG""4230,
PROJECTIONTransverse_Mercator,
PARAMETERlatitude_of_origin,
PARAMETERcentral_meridian,
PARAMETERscale_factor,
PARAMETERfalse_easting,
PARAMETERfalse_northing,
UNIT["metre",1,
AUTHORITYEPSG""9001,
AUTHORITYEPSG""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

History

#1 Updated by mmassing - over 10 years ago

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

#2 Updated by mmassing - over 10 years ago

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.

#3 Updated by alobo - over 10 years ago

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

#4 Updated by mmassing - over 10 years ago

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

#5 Updated by rbivand - over 10 years ago

Replying to [comment:2 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.

#6 Updated by rbivand - over 10 years ago

Replying to [comment:5 rbivand]:

Replying to [comment:2 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.

#7 Updated by Marco Hugentobler over 10 years ago

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?

#8 Updated by rbivand - over 10 years ago

Replying to [comment:7 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).

#9 Updated by alobo - over 10 years ago

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

#10 Updated by alobo - over 10 years ago

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

#11 Updated by Giuseppe Sucameli over 10 years ago

I think there is not only one issue.

First issue:

Replying to [comment:10 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 [comment:10 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)?

#12 Updated by alobo - over 10 years ago

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

#13 Updated by Giovanni Manghi almost 9 years ago

  • Target version changed from Version 1.7.0 to Version 1.7.4

#14 Updated by Paolo Cavallini over 8 years ago

  • Affected QGIS version set to master
  • Crashes QGIS or corrupts data set to No
  • Target version changed from Version 1.7.4 to Version 1.8.0

#15 Updated by Paolo Cavallini about 8 years ago

  • Target version changed from Version 1.8.0 to Version 2.0.0

#16 Updated by Giovanni Manghi almost 8 years ago

  • Status changed from Open to Feedback
  • Assignee deleted (nobody -)
  • Status info deleted (1)
  • Pull Request or Patch supplied set to No

Can sample data be provided again? The links are not working anymore.

#17 Updated by alobo - almost 8 years ago

Updated link:
http://dl.dropbox.com/u/3180464/test2.tif

QGIS raster/properties/meteadata is now:
GDAL provider
VRT
Virtual Raster
Dataset Description
Band 1
Dimensions:
X: 316 Y: 177 Bands: 1
Origin:
424389,4.68549e+06
Pixel Size:
279.887,-279.887
No Data Value
2.14748e+09
Data Type:
GDT_Float64 - Sixty four bit floating point

gdalinfo has not changed (note differences in pixel size)

Agus

#18 Updated by Giovanni Manghi almost 8 years ago

  • Operating System deleted (Debian)
  • Status changed from Feedback to Open

#19 Updated by Jürgen Fischer about 6 years ago

  • Target version changed from Version 2.0.0 to Future Release - Lower Priority

#20 Updated by Giovanni Manghi over 3 years ago

  • Regression? set to No
  • Easy fix? set to No

#21 Updated by Giovanni Manghi over 1 year ago

  • Status changed from Open to Closed
  • Resolution set to end of life

Also available in: Atom PDF