How to calculate affine transformation by Python gdal as same as gdal_translate?

Question

I've made a python script which tries to reproduce what gdal_translate does (given a bounding box BBOX in EPSG:4326 and an UTM32 .jp2 file. Then, translates the BBOX from EPSG:4326 to EPSG:32632 (UTM32 North) and subsets the result from the geotiff file.

Gdalinfo of the Sentinel 2 UTM32 file:

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)

Below is the python script:

from osgeo import gdal
import sys
import ogr, osr
import math
import os

utm32_file = "T32UQD_20200421T102021_TCI_60m.jp2"
gdalwarp_crop_file = "gdal_cropped.tiff"

# GIVEN BBOX in EPSG:4326
X_MIN, Y_MIN, X_MAX, Y_MAX = 12.46826, 52.58852, 12.48516, 52.60034

def project(pointX, pointY):
    # Spatial Reference System
    inputEPSG = 4326
    outputEPSG = 32632

    # create a geometry from coordinates
    point = ogr.Geometry(ogr.wkbPoint)
    point.AddPoint(pointX, pointY)

    # create coordinate transformation
    inSpatialRef = osr.SpatialReference()
    inSpatialRef.ImportFromEPSG(inputEPSG)

    outSpatialRef = osr.SpatialReference()
    outSpatialRef.ImportFromEPSG(outputEPSG)

    coordTransform = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)

    # transform point
    point.Transform(coordTransform)

    # print point in EPSG 4326
    return point.GetX(), point.GetY()

xmin, ymin = project(X_MIN, Y_MIN)
xmax, ymax = project(X_MAX, Y_MAX)

# coordinates to get pixel values for (as tuples of points)
points = [(xmin, ymin), (xmax, ymax)]

# open the raster file
ds = gdal.Open(utm32_file)

if ds is None:
    print('Could not open the raster file')
    sys.exit(1)
else:
    print('The raster file was opened satisfactorily')

# get georeference info
transform = ds.GetGeoTransform() # (-2493045.0, 30.0, 0.0, 3310005.0, 0.0, -30.0)
xOrigin = transform[0] # -2493045.0
yOrigin = transform[3] # 3310005.0
pixelWidth = transform[1] # 30.0
pixelHeight = transform[5] # -30.0

band = ds.GetRasterBand(1) # 1-based index

data = band.ReadAsArray()

grid_x = 0
grid_y = 0

# loop through the coordinates
for point in points:
    x = point[0]
    y = point[1]

    xOffset = int((x - xOrigin) / pixelWidth)
    yOffset = int((y - yOrigin) / pixelHeight)

    if grid_x == 0:
        grid_x = xOffset
        grid_y = yOffset
    else:
        grid_x = abs(xOffset - grid_x) + 1
        grid_y = abs(yOffset - grid_y) + 1

gdal_subset = str(points[0][0]) + " " + str(points[1][1]) + " " + str(points[1][0]) + " " + str(points[0][1])
command = "gdal_translate -projwin " + gdal_subset + " " + utm32_file + " " + gdalwarp_crop_file + " >> nohup.out"
os.system(command)

# PRINT GRID SIZES (width x height)
print "python:    ", grid_x, grid_y

ds = gdal.Open(gdalwarp_crop_file)
print "gdal:      ", ds.RasterXSize, ds.RasterYSize

The problem is the affine transformation in the script sometimes returns correct grid sizes (width x height) as gdal_translate, but in many cases it returns more few pixels than gdal_translate. For example output of the python script is:

python:     19 x 24
gdal:       18 x 23

What could be the problem here to return the same result (grid sizes) as gdal_translate?

Answer 1

Your original question was how to calculate the same number of pixels as gdal_translate. I honestly do not know the answer to that question other than by looking at the gdal_translate code, where there are several nuances in cropping the source raster. You could compare the gdal_translate code with your code and possibly step through both codes side by side in a debugger.

I guess my approach to solve this problem would not be to attempt to rewrite gdal_translate, but rather to use the gdal.Translate python method as described in the answer to this question.