Difference results from QGIS project on fly and gdalwarp for EPSG:4326 of Sentinel 2 UTM zone 32?

Question

I always thought QGIS and gdalwarp should return same result for projecting a 2D Sentinel 2 image (e.g: T32UQD_20200421T102021_TCI_60m.jp2) from UTM 32 to EPSG:4326.

What I did was:

• Open the .jp2 in QGis 2.14.9 with 'Enable on the fly CRS transformation'.

• Warp the .jp2 to geotiff by gdalwarp:

  gdalwarp T32UQD_20200421T102021_TCI_60m.jp2 -t_srs EPSG:4326 4.warped.tiff

But the output results from them look different significantly (on the left hand side with QGIS it looks rotated and on the right hand side with gdalwarp it looks aligned).

So, what should be the correct result?

gdalinfo T32UQD_20200421T102021_TCI_60m.jp2

Size is 1830, 1830
Coordinate System is:
PROJCS["WGS 84 / UTM zone 32N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AXIS["Latitude",NORTH],
        AXIS["Longitude",EAST],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",9],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","32632"]]
Origin = (699960.000000000000000,5900040.000000000000000)
Pixel Size = (60.000000000000000,-60.000000000000000)
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  699960.000, 5900040.000) ( 11d59'40.71"E, 53d12'43.17"N)
Lower Left  (  699960.000, 5790240.000) ( 11d55'40.15"E, 52d13'34.41"N)
Upper Right (  809760.000, 5900040.000) ( 13d38' 3.06"E, 53d 9'33.92"N)
Lower Right (  809760.000, 5790240.000) ( 13d31'51.72"E, 52d10'31.74"N)
Center      (  754860.000, 5845140.000) ( 12d46'18.91"E, 52d41'45.89"N)

and warped by gdalwarp: 4.warped.tiff

Size is 2383, 1447
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (11.927820031375013,53.211992871316582)
Pixel Size = (0.000716123003772,-0.000716123003772)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  11.9278200,  53.2119929) ( 11d55'40.15"E, 53d12'43.17"N)
Lower Left  (  11.9278200,  52.1757629) ( 11d55'40.15"E, 52d10'32.75"N)
Upper Right (  13.6343411,  53.2119929) ( 13d38' 3.63"E, 53d12'43.17"N)
Lower Right (  13.6343411,  52.1757629) ( 13d38' 3.63"E, 52d10'32.75"N)
Center      (  12.7810806,  52.6938779) ( 12d46'51.89"E, 52d41'37.96"N)

Answer 1

I don't see anything wrong here. Gdalwarp creates a new raster with cellsize in units of the target CRS sampling from the original raster values as needed, while OTF reprojection just squeezes the original raster to the new CRS without recalculating cell values.

You see the difference in the corner coordinates. The gdalwarp output is a rectangle in the target CRS, while the source raster is a rectangle in the source CRS.

Answer 2

To elaborate on @AndreJ's answer, which correctly said that this is entirely expected: It is also a good illustration of what "on the fly" means. The original UTM32N raster remains the same, the pixel corner locations remain the same. It's just that their coordinates are on the fly translated to WGS84.

As you can see from the gdalinfo output, the two rasters have different numbers of rows and columns. The UTM32N one has square pixels in UTM32N, measured in meters (60x60m), and each pixel has the same area. These will look like distorted rectangles when represented in WGS84. Conversely, the reprojected (warped) raster in WGS84 has square pixels in degrees. These look square in WGS84, but wouldn't be square when measured in meters. Their area varies.

The warping was done so as to roughly preserve, I think, pixel area. The visual differences between the raster values are from resampling.

QGIS and gdalwarp should produce the same result if you ask them to do the same thing (ie, produce the same warped raster - not just the same target CRS, but the same origin, pixel size, and resampling method. But not when they are doing different things.