QGIS » QGIS Application

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

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.

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

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

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.

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.

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?

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

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

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

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]]

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

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

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