How to prevent gdalwarp from swapping axes?

Question

Executing gdalwarp on a file which has an axis mapping (Data axis to CRS axis mapping: 2,1) results in a transposed image (axes are swapped).

Is there a way to prevent this? I tried to use t_srs option, to no avail.

I know I could use, for example, Python bindings and transpose the data with numpy. But I'd rather having a solution with GDAL's CLI.

The command I use:

gdalwarp -of GTiff input.JP2 output.tif

See the sizes in gdalinfo outputs.

JPEG2000 input file's gdalinfo:

Size is 2407, 3729
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
GeoTransform =
  4.424925925925925, 0, 3.7037037037e-05
  50.9246111111111, -1.8518518519e-05, 0
Metadata:
  METADATATYPE=DIMAP
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (   4.4249259,  50.9246111) (  4d25'29.73"E, 50d55'28.60"N)
Lower Left  (   4.5630370,  50.9246111) (  4d33'46.93"E, 50d55'28.60"N)
Upper Right (   4.4249259,  50.8800370) (  4d25'29.73"E, 50d52'48.13"N)
Lower Right (   4.5630370,  50.8800370) (  4d33'46.93"E, 50d52'48.13"N)
Center      (   4.4939815,  50.9023241) (  4d29'38.33"E, 50d54' 8.37"N)
Band 1 Block=2048x2048 Type=UInt16, ColorInterp=Red
  Overviews: 1204x1865, 602x933, 301x467, 151x234, 76x117
  Overviews: arbitrary
  Image Structure Metadata:
    COMPRESSION=JPEG2000
    NBITS=12
Band 2 Block=2048x2048 Type=UInt16, ColorInterp=Green
  Overviews: 1204x1865, 602x933, 301x467, 151x234, 76x117
  Overviews: arbitrary
  Image Structure Metadata:
    COMPRESSION=JPEG2000
    NBITS=12
Band 3 Block=2048x2048 Type=UInt16, ColorInterp=Blue
  Overviews: 1204x1865, 602x933, 301x467, 151x234, 76x117
  Overviews: arbitrary
  Image Structure Metadata:
    COMPRESSION=JPEG2000
    NBITS=12
Band 4 Block=2048x2048 Type=UInt16, ColorInterp=Undefined
  Overviews: 1204x1865, 602x933, 301x467, 151x234, 76x117
  Overviews: arbitrary
  Image Structure Metadata:
    COMPRESSION=JPEG2000
    NBITS=12

GeoTIFF output's gdalinfo:

Size is 4224, 1363
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (4.424925925925925,50.924611111111105)
Pixel Size = (0.000032698049453,-0.000032698049453)
Metadata:
  AREA_OR_POINT=Area
  METADATATYPE=DIMAP
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (   4.4249259,  50.9246111) (  4d25'29.73"E, 50d55'28.60"N)
Lower Left  (   4.4249259,  50.8800437) (  4d25'29.73"E, 50d52'48.16"N)
Upper Right (   4.5630425,  50.9246111) (  4d33'46.95"E, 50d55'28.60"N)
Lower Right (   4.5630425,  50.8800437) (  4d33'46.95"E, 50d52'48.16"N)
Center      (   4.4939842,  50.9023274) (  4d29'38.34"E, 50d54' 8.38"N)
Band 1 Block=4224x1 Type=UInt16, ColorInterp=Red
Band 2 Block=4224x1 Type=UInt16, ColorInterp=Green
Band 3 Block=4224x1 Type=UInt16, ColorInterp=Blue
Band 4 Block=4224x1 Type=UInt16, ColorInterp=Undefined

Answer 1

The georeferencing of the source image is not most typical. Look at the corner coordinates first

Corner Coordinates:
Upper Left  (   4.4249259,  50.9246111) (  4d25'29.73"E, 50d55'28.60"N)
Lower Left  (   4.5630370,  50.9246111) (  4d33'46.93"E, 50d55'28.60"N)
Upper Right (   4.4249259,  50.8800370) (  4d25'29.73"E, 50d52'48.13"N)
Lower Right (   4.5630370,  50.8800370) (  4d33'46.93"E, 50d52'48.13"N)

For non-rotated images the longitude of upper left and lower left are the same but for this image they are 4.4249259 and 4.5630370.

Then have a look at the GeoTransform

GeoTransform =
  4.424925925925925, 0, 3.7037037037e-05
  50.9246111111111, -1.8518518519e-05, 0

Compare the parameters with the documentation https://gdal.org/user/raster_data_model.html

Index of the parameters is

0,1,2
3,4,5

Parameters 0 and 3 point the origin and they are fine. Normally parameters 1 and 5 show the pixel size in west-east and north-south directions. In your images they are both zero. Parameters 2 and 4 are zero for north-up images but you have values 3.7037037037e-05 and -1.8518518519e-05.

This means that your original is an east-up (or west-up) image. Gdalwarp is producing north-up images and the result that you get is expected.

Gdal_translate can change the format but it does not warp and therefore it probably reserves the uncommon georeferencing. Have a try with

gdal_translate -of GTiff input.jp2 output.tif

If you think that the original image actually is north-up image then contact your image provider and ask them to correct the metadata of their image. Without having access to your JPEG2000 image I can't say what is wrong but I guess that the GML coverage that is stored inside the JPEG2000 file is using wrong axis order. If you can't make image provider to fix their data then then there could be some workaround by using GDAL. That is worth another question, though.