search

qissue-bot / QGIS

QGIS is a free, open source, cross platform (lin/win/mac) geographical information system (GIS)
http://qgis.org
GNU General Public License v2.0
0 stars 0 forks source link

Discrepancy between gdalinfo and QGIS raster metadata #2431

Open qissue-bot opened 5 years ago

qissue-bot commented 5 years ago

Author Name: alobo - (alobo -) Original Redmine Issue: 2444, https://issues.qgis.org/issues/2444

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]], UNIT[[degree]], AUTHORITY[[EPSG""4230]], PROJECTION[[Transverse_Mercator]], PARAMETER[[latitude_of_origin]], PARAMETER[[central_meridian]], PARAMETER[[scale_factor]], PARAMETER[[false_easting]], PARAMETER[[false_northing]], 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

qissue-bot commented 5 years ago

Original Redmine Comment Author Name: mmassing - (mmassing -) Original Date: 2010-02-15T03:03:14.000Z

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

qissue-bot commented 5 years ago

Original Redmine Comment Author Name: mmassing - (mmassing -) Original Date: 2010-02-15T03:13:11.000Z

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.

qissue-bot commented 5 years ago

Original Redmine Comment Author Name: alobo - (alobo -) Original Date: 2010-02-15T04:50:05.000Z

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

qissue-bot commented 5 years ago

Original Redmine Comment Author Name: mmassing - (mmassing -) Original Date: 2010-02-15T05:06:17.000Z

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

qissue-bot commented 5 years ago

Original Redmine Comment Author Name: rbivand - (rbivand -) Original Date: 2010-02-15T05:20:53.000Z

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.

qissue-bot commented 5 years ago

Original Redmine Comment Author Name: rbivand - (rbivand -) Original Date: 2010-02-15T07:07:30.000Z

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.

qissue-bot commented 5 years ago

Original Redmine Comment Author Name: Marco Hugentobler (Marco Hugentobler) Original Date: 2010-02-17T02:01:52.000Z

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?

qissue-bot commented 5 years ago

Original Redmine Comment Author Name: rbivand - (rbivand -) Original Date: 2010-02-17T02:29:56.000Z

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).

qissue-bot commented 5 years ago

Original Redmine Comment Author Name: alobo - (alobo -) Original Date: 2010-02-17T02:35:10.000Z

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

qissue-bot commented 5 years ago

Original Redmine Comment Author Name: alobo - (alobo -) Original Date: 2010-02-17T07:45:25.000Z

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

qissue-bot commented 5 years ago

Original Redmine Comment Author Name: Giuseppe Sucameli (Giuseppe Sucameli) Original Date: 2010-05-03T13:18:44.000Z

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:

  1. 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)?
qissue-bot commented 5 years ago

Original Redmine Comment Author Name: alobo - (alobo -) Original Date: 2010-05-03T21:04:23.000Z

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

qissue-bot commented 5 years ago

Original Redmine Comment Author Name: Giovanni Manghi (Giovanni Manghi) Original Date: 2012-10-05T05:11:46.000Z

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

qissue-bot commented 5 years ago

Original Redmine Comment Author Name: alobo - (alobo -) Original Date: 2012-10-06T20:53:56.000Z

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