Crop GTiff to match bounding box

Question

Related to a previous question, I need to crop my dataset for visualization purposes. By crop I mean getting rid of the nodata values and rotating the result.

My raster originally looks like this:

I am able to export the raster as a PNG using gdal_translate as in gdal_translate -of PNG downloads\2023-05-01.tiff test.png -scale -ot Byte but this returns the image as a rotated PNG.

I would like to undo this geometric transformation and obtain something like this even if it means losing information, georreferentiation or resampling:

Not sure if this can be done only with GDAL, but I'm using the Python bindings, so I have much more tools available besides GDAL if not possible. I'm unsure how this operation is called, so I don't know where or what to look for in the docs.

Running gdalinfo for my dataset as gdalinfo -stats downloads/2023-05-01.tiff outputs the following:

Driver: GTiff/GeoTIFF
Files: downloads/2023-05-01.tiff
Size is 1818, 1593
Coordinate System is:
PROJCRS["WGS 84 / UTM zone 34S",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 34S",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",21,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",10000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Navigation and medium accuracy spatial referencing."],
        AREA["Between 18┬░E and 24┬░E, southern hemisphere between 80┬░S and equator, onshore and offshore. Angola. Botswana. Democratic Republic of the Congo (Zaire). Namibia. South Africa. Zambia."],
        BBOX[-80,18,0,24]],
    ID["EPSG",32734]]
Data axis to CRS axis mapping: 1,2
Origin = (755100.000000000000000,7811940.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  LAYOUT=COG
  COMPRESSION=LZW
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  755100.000, 7811940.000) ( 23d26' 4.68"E, 19d46'18.38"S)
Lower Left  (  755100.000, 7796010.000) ( 23d26'12.58"E, 19d54'56.15"S)
Upper Right (  773280.000, 7811940.000) ( 23d36'28.84"E, 19d46' 9.57"S)
Lower Right (  773280.000, 7796010.000) ( 23d36'37.30"E, 19d54'47.28"S)
Center      (  764190.000, 7803975.000) ( 23d31'20.85"E, 19d50'32.92"S)
Band 1 Block=512x512 Type=UInt16, ColorInterp=Gray
  Minimum=1077.000, Maximum=13880.000, Mean=2018.113, StdDev=351.323
  NoData Value=0
  Overviews: 909x796, 454x398, 227x199
  Offset: -0.1,   Scale:0.0001
  Metadata:
    STATISTICS_MINIMUM=1077
    STATISTICS_MAXIMUM=13880
    STATISTICS_MEAN=2018.1125217712
    STATISTICS_STDDEV=351.32337310468
    STATISTICS_VALID_PERCENT=55.21
Band 2 Block=512x512 Type=UInt16, ColorInterp=Undefined
  Minimum=1073.000, Maximum=12488.000, Mean=1811.102, StdDev=233.957
  NoData Value=0
  Overviews: 909x796, 454x398, 227x199
  Offset: -0.1,   Scale:0.0001
  Metadata:
    STATISTICS_MINIMUM=1073
    STATISTICS_MAXIMUM=12488
    STATISTICS_MEAN=1811.1022441941
    STATISTICS_STDDEV=233.95662853744
    STATISTICS_VALID_PERCENT=55.21
Band 3 Block=512x512 Type=UInt16, ColorInterp=Undefined
  Minimum=1065.000, Maximum=9832.000, Mean=1572.875, StdDev=191.223
  NoData Value=0
  Overviews: 909x796, 454x398, 227x199
  Offset: -0.1,   Scale:0.0001
  Metadata:
    STATISTICS_MINIMUM=1065
    STATISTICS_MAXIMUM=9832
    STATISTICS_MEAN=1572.8751681504
    STATISTICS_STDDEV=191.22301182326
    STATISTICS_VALID_PERCENT=55.21

Answer 1

Saving into PNG format does not actually return a roteted image. Even the original TIFF image is a rectangle with the image area that appears rotated at the centre and nodata triangles at the corners. The image file is georeferenced including the nodata areas. If you look at the Corner Coordinates in the gdalinfo report you will notice that the coordinates present the corners of the full image file, not the corners of the data area. Nodata aware viewer just knows to make the nodata pixels transparent.

Therefore it is not possible to "undo this geometric transformation" as you suggest. What you can do is to select the area that covers the real image data (often called as "footprint"), crop the image, and stretch the image from the corners.

The GDAL workflow could be:

• Find the four corners and insert ground control points (GCP) into the TIFF for example with gdal_translate https://gdal.org/en/latest/programs/gdal_translate.html#cmdoption-gdal_translate-gcp. Use some fake coordinates that makes a rectangle as easting and northing (for example (0,0), (0,1000), (1000,1000), and (1000,0))
• Warp and crop the image with gdalwarp. Use "target extents" for cropping https://gdal.org/en/latest/programs/gdalwarp.html#cmdoption-gdalwarp-te.

Another option would be to use some generic image processing software like PhotoShop or GIMP https://graphicdesign.stackexchange.com/questions/137304/how-to-fix-photo-perspective-and-crop-to-a-perfect-rectangle. That would be easy to do manually but hard to automate. For the automatic processing you should anyway find the corners of the image footprint automatically. It is possible that the satellite image producer has included the footprint in some metadata document.