Opening raster file ordered bottom-up with non-square pixels produces wrong output

[jjimenezshaw]

What is the bug or the crash?

Opening a raster file with lines ordered bottom-up, with a non-squre pixel size produces something wrong. The pixel sizes shown in the properties is something in between the two, and the image displayed in the screen is wrong.
The same data ordered top-down is correctly managed.

This kind of non square pixel rasters can be normal on geographic coordinates. That way, they represent something "more square" in the field.

Steps to reproduce the issue

I created this XYZ file manually (bottom-up.xyz)

y x z
30.0 10.0 1
30.0 10.2 2
30.0 10.4 3
30.1 10.0 4
30.1 10.2 5
30.1 10.4 6
30.2 10.0 7
30.2 10.2 8
30.2 10.4 9

with 0.1 step in Y (latitude) and 0.2 in X (longitude)

$ gdalinfo bottom-up.xyz
Driver: XYZ/ASCII Gridded XYZ
Files: bottom-up.xyz
Size is 3, 3
Origin = (9.900000000000000,29.949999999999999)
Pixel Size = (0.200000000000000,0.100000000000000)
Corner Coordinates:
Upper Left  (   9.9000000,  29.9500000) 
Lower Left  (   9.9000000,  30.2500000) 
Upper Right (  10.5000000,  29.9500000) 
Lower Right (  10.5000000,  30.2500000) 
Center      (  10.2000000,  30.1000000) 
Band 1 Block=3x1 Type=Byte, ColorInterp=Undefined
  Min=1.000 Max=9.000

The pixel size displayed in the properties in QGIS is 0.1581138830084189706,-0.1581138830084189706.

I can easily reorder the values to get a top-down file with the command

cat bottom-up.xyz | sort -s -k1,1r -k2,2  > top-down.xyz

producing this file

y x z
30.2 10.0 7
30.2 10.2 8
30.2 10.4 9
30.1 10.0 4
30.1 10.2 5
30.1 10.4 6
30.0 10.0 1
30.0 10.2 2
30.0 10.4 3

$ gdalinfo top-down.xyz 
Driver: XYZ/ASCII Gridded XYZ
Files: top-down.xyz
Size is 3, 3
Origin = (9.900000000000000,30.250000000000000)
Pixel Size = (0.200000000000000,-0.100000000000000)
Corner Coordinates:
Upper Left  (   9.9000000,  30.2500000) 
Lower Left  (   9.9000000,  29.9500000) 
Upper Right (  10.5000000,  30.2500000) 
Lower Right (  10.5000000,  29.9500000) 
Center      (  10.2000000,  30.1000000) 
Band 1 Block=3x1 Type=Byte, ColorInterp=Undefined
  Min=1.000 Max=9.000

The pixel size displayed in the properties in QGIS is 0.2000000000000001776,-0.09999999999999964473

This file top-down.xyz is displayed properly in QGIS.

I tried ordering the columns as lon-lat-z, with the same results.

With same spacing in both directions no problem is detected in the display.

Versions

(I tested also with conda version 3.34.1-Prizren, with GDAL 3.8.0)

<style type="text/css"> p, li { white-space: pre-wrap; } </style>

QGIS version	3.34.3-Prizren	QGIS code revision	4737323
Qt version	5.15.3
Python version	3.10.12
GDAL/OGR version	3.4.1
PROJ version	8.2.1
EPSG Registry database version	v10.041 (2021-12-03)
GEOS version	3.10.2-CAPI-1.16.0
SQLite version	3.37.2
PDAL version	2.3.0
PostgreSQL client version	14.10 (Ubuntu 14.10-0ubuntu0.22.04.1)
SpatiaLite version	5.0.1
QWT version	6.1.4
QScintilla2 version	2.11.6
OS version	Ubuntu 22.04.3 LTS
Active Python plugins
FreehandRasterGeoreferencer	0.8.3
db_manager	0.1.20
processing	2.12.99
grassprovider	2.12.99
MetaSearch	0.3.6

QGIS version 3.34.3-Prizren QGIS code revision [4737323](https://github.com/qgis/QGIS/commit/47373234acd) Qt version 5.15.3 Python version 3.10.12 GDAL/OGR version 3.4.1 PROJ version 8.2.1 EPSG Registry database version v10.041 (2021-12-03) GEOS version 3.10.2-CAPI-1.16.0 SQLite version 3.37.2 PDAL version 2.3.0 PostgreSQL client version 14.10 (Ubuntu 14.10-0ubuntu0.22.04.1) SpatiaLite version 5.0.1 QWT version 6.1.4 QScintilla2 version 2.11.6 OS version Ubuntu 22.04.3 LTS

Active Python plugins
FreehandRasterGeoreferencer
0.8.3
db_manager
0.1.20
processing
2.12.99
grassprovider
2.12.99
MetaSearch
0.3.6

Supported QGIS version

• I'm running a supported QGIS version according to the roadmap.

New profile

• I tried with a new QGIS profile

Additional context

Based on the difference between gdalinfo and the properties, I think the problem could be in QGIS. In case you think it is a problem in GDAL, I can open a ticket there.

[agiudiceandrea]

Moreover, the bottom-up.xyz layer has the incorrect property of "Dimensions X: 4 Y: 2" in QGIS.

[jjimenezshaw]

I am trying to debug the code. But maybe the reason is in

QGIS/src/core/providers/gdal/qgsgdalprovider.cpp

Line 3590 in 3a18496

// Check if we need a warped VRT for this file.

  // Check if we need a warped VRT for this file.
  bool hasGeoTransform = GDALGetGeoTransform( mGdalBaseDataset, mGeoTransform ) == CE_None;
  if ( ( hasGeoTransform
         && ( mGeoTransform[1] < 0.0
              || mGeoTransform[2] != 0.0
              || mGeoTransform[4] != 0.0
              || mGeoTransform[5] > 0.0 ) )
       || GDALGetGCPCount( mGdalBaseDataset ) > 0
       || GDALGetMetadata( mGdalBaseDataset, "RPC" ) )
  {
    QgsDebugMsgLevel( QStringLiteral( "Creating Warped VRT." ), 2 );

If that is the case (I am not sure), the trigger would be mGeoTransform[5] > 0.0. But then the VRT is setting an square pixel. That distorts the final output. In that case, at least I would expect a message telling the user that some data is not true.

[elpaso]

Working on a GDAL fix.