How to use the gdal.ApplyGeoTransform function in GDAL

# open the GDAL file
sourceImage = gdal.Open(imageFileName, gdal.GA_ReadOnly)
rpcMetaData = sourceImage.GetMetadata('RPC')
myarray = np.array(sourceImage.GetRasterBand(1).ReadAsArray())
myarray = myarray*0+255

gt = sourceImage.GetGeoTransform() # captures origin and pixel size
print('Origin:', (gt[0], gt[3]))
print('Pixel size:', (gt[1], gt[5]))

gdal.ApplyGeoTransform(gt,0,0)

left = gdal.ApplyGeoTransform(gt,0,0)[0]
top = gdal.ApplyGeoTransform(gt,0,0)[1]
right = gdal.ApplyGeoTransform(gt,sourceImage.RasterXSize,sourceImage.RasterYSize)[0]
bottom = gdal.ApplyGeoTransform(gt,sourceImage.RasterXSize,sourceImage.RasterYSize)[1]

print(left, top, right, bottom)
osmFileName = '../data/map.osm'
#osmFileName = '../data/state/cb_2016_us_state_20m.shp'

ds = ogr.Open(osmFileName)
layer = ds.GetLayerByIndex(3)

nameList = []
for feature in layer:
    if feature.GetField("tourism") != None or feature.GetField("building") != None:  # only streets
        flag = True
        #print("haha")
        name = feature.GetField("name")
        geom = feature.GetGeometryRef()

#----------------------------

        # check raster parameters to find default zoom range
        ld('automatic zoom levels')

        # modify target srs to allow charts crossing meridian 180
        shifted_srs = self.shift_srs()

        t_ds = gdal.AutoCreateWarpedVRT(self.src_ds, None, txt2wkt(shifted_srs))
        geotr = t_ds.GetGeoTransform()
        res = (geotr[1], geotr[5])
        max_zoom = max(self.res2zoom_xy(res))

        # calculate min_zoom
        ul_c = (geotr[0], geotr[3])
        lr_c = gdal.ApplyGeoTransform(geotr, t_ds.RasterXSize, t_ds.RasterYSize)
        wh = (lr_c[0]-ul_c[0], ul_c[1]-lr_c[1])
        ld('ul_c, lr_c, wh', ul_c, lr_c, wh)
        min_zoom = min(self.res2zoom_xy([wh[i]/abs(self.tile_dim[i]) for i in (0, 1)]))

        self.set_zoom_range(zoom_parm, (min_zoom, max_zoom))

scale = args.scale

    inputImage = gdal_utils.gdal_open(args.input_image, gdal.GA_ReadOnly)
    band = inputImage.GetRasterBand(1)
    if (not band.DataType == gdal.GDT_Byte):
        raise RuntimeError(
            "Input image {} does not have Byte type. Use msi-to-rgb.py to-8bit.py "
            "to convert it.".format(args.input_image))

    projection = inputImage.GetProjection()
    inputImageSrs = osr.SpatialReference(projection)
    gt = inputImage.GetGeoTransform()  # captures origin and pixel size

    left, top = gdal.ApplyGeoTransform(gt, 0, 0)
    right, bottom = gdal.ApplyGeoTransform(gt, inputImage.RasterXSize, inputImage.RasterYSize)
    band = None

    print("Resize and edge detection: {}".format(os.path.basename(args.input_image)))
    color_image = cv2.imread(args.input_image)
    small_color_image = numpy.zeros(
        (int(color_image.shape[0]*scale),
         int(color_image.shape[1]*scale), 3), dtype=numpy.uint8)
    if scale != 1.0:
        small_color_image = cv2.resize(color_image, None, fx=scale, fy=scale)
        color_image = small_color_image
    grayimg = cv2.cvtColor(color_image, cv2.COLOR_BGR2GRAY)
    edge_img = cv2.Canny(grayimg, 100, 200)
    if args.debug:
        cv2.imwrite(os.path.splitext(args.output_mask)[0] + '_edge.tif', edge_img)

    model = {}

import gdal, ogr, os, osr
import numpy as np


ds = gdal.Open('../data/Dayton.ntf')
print('File list:', ds.GetFileList())
print('Width:', ds.RasterXSize)
print('Height:', ds.RasterYSize)
print('Coordinate system:', ds.GetProjection())
gt = ds.GetGeoTransform() # captures origin and pixel size
print('Origin:', (gt[0], gt[3]))
print('Pixel size:', (gt[1], gt[5]))
print('Upper Left Corner:', gdal.ApplyGeoTransform(gt,0,0))
print('Upper Right Corner:', gdal.ApplyGeoTransform(gt,ds.RasterXSize,0))
print('Lower Left Corner:', gdal.ApplyGeoTransform(gt,0,ds.RasterYSize))
print('Lower Right Corner:',
gdal.ApplyGeoTransform(gt,ds.RasterXSize,ds.RasterYSize))
print('Center:', gdal.ApplyGeoTransform(gt,ds.RasterXSize/2,ds.RasterYSize/2))
print('Metadata:', ds.GetMetadata())
print('Image Structure Metadata:', ds.GetMetadata('IMAGE_STRUCTURE'))
print('Number of bands:', ds.RasterCount)

for i in range(1, ds.RasterCount+1):
    band = ds.GetRasterBand(i) # in GDAL, band are indexed starting at 1!
    interp = band.GetColorInterpretation()
    interp_name = gdal.GetColorInterpretationName(interp)
    (w,h)=band.GetBlockSize()
    print('Band %d, block size %dx%d, color interp %s' % (i,w,h,interp_name))

#----------------------------

        # check raster parameters to find default zoom range
        ld('automatic zoom levels')

        # modify target srs to allow charts crossing meridian 180
        shifted_srs = self.shift_srs()

        t_ds = gdal.AutoCreateWarpedVRT(self.src_ds, None, txt2wkt(shifted_srs))
        geotr = t_ds.GetGeoTransform()
        res = (geotr[1], geotr[5])
        max_zoom = max(self.res2zoom_xy(res))

        # calculate min_zoom
        ul_c = (geotr[0], geotr[3])
        lr_c = gdal.ApplyGeoTransform(geotr, t_ds.RasterXSize, t_ds.RasterYSize)
        wh = (lr_c[0]-ul_c[0], ul_c[1]-lr_c[1])
        ld('ul_c, lr_c, wh', ul_c, lr_c, wh)
        min_zoom = min(self.res2zoom_xy([wh[i]/abs(self.tile_dim[i]) for i in (0, 1)]))

        self.set_zoom_range(zoom_parm, (min_zoom, max_zoom))

:param lat: latitude
    :param lon: longitude
    :param transform: GDAL GeoTransform
    :return: x, y pair
    """
    # Get the reverse direction GeoTransform
    _, transform = gdal.InvGeoTransform(transform)
    
    '''
    ulx, uly = transform[0], transform[3]
    psx, psy = transform[1], transform[5]
    x = int(math.floor(ulx + psx * lon))
    y = int(math.floor(uly + psy * lat))
    '''
    # This is to fix to make nearest rounding before casting to integer and apply gdal utility
    x, y= gdal.ApplyGeoTransform(transform, lon, lat)
    return int(round(x)), int(round(y))

#lowThresh = 0.5*high_thresh
#edge_img = cv2.Canny(grayimg, lowThresh, high_thresh)
#edge_img = cv2.Canny(grayimg, 200, 300)
edge_img = cv2.Canny(grayimg, 100, 200)
cv2.imwrite('../data/{}_edge.jpg'.format(basename),edge_img)

# open the GDAL file
sourceImage = gdal.Open(input_img, gdal.GA_ReadOnly)
rpcMetaData = sourceImage.GetMetadata('RPC')
full_res_mask = np.zeros((sourceImage.RasterYSize,sourceImage.RasterXSize,4))

gt = sourceImage.GetGeoTransform() # captures origin and pixel size
print('Origin:', (gt[0], gt[3]))
print('Pixel size:', (gt[1], gt[5]))

left = gdal.ApplyGeoTransform(gt,0,0)[0]
top = gdal.ApplyGeoTransform(gt,0,0)[1]
right = gdal.ApplyGeoTransform(gt,sourceImage.RasterXSize,sourceImage.RasterYSize)[0]
bottom = gdal.ApplyGeoTransform(gt,sourceImage.RasterXSize,sourceImage.RasterYSize)[1]

model = {}
model['corners'] = [left, top, right, bottom]
model['project_model'] = gt 
model['scale'] = scale

ds = ogr.Open(args.input_osm)
if ds.GetLayerCount()<=1:
    layer = ds.GetLayerByIndex(0)
else:
    layer = ds.GetLayerByIndex(3)
nameList = []

import numpy as np


ds = gdal.Open('../data/Dayton.ntf')
print('File list:', ds.GetFileList())
print('Width:', ds.RasterXSize)
print('Height:', ds.RasterYSize)
print('Coordinate system:', ds.GetProjection())
gt = ds.GetGeoTransform() # captures origin and pixel size
print('Origin:', (gt[0], gt[3]))
print('Pixel size:', (gt[1], gt[5]))
print('Upper Left Corner:', gdal.ApplyGeoTransform(gt,0,0))
print('Upper Right Corner:', gdal.ApplyGeoTransform(gt,ds.RasterXSize,0))
print('Lower Left Corner:', gdal.ApplyGeoTransform(gt,0,ds.RasterYSize))
print('Lower Right Corner:',
gdal.ApplyGeoTransform(gt,ds.RasterXSize,ds.RasterYSize))
print('Center:', gdal.ApplyGeoTransform(gt,ds.RasterXSize/2,ds.RasterYSize/2))
print('Metadata:', ds.GetMetadata())
print('Image Structure Metadata:', ds.GetMetadata('IMAGE_STRUCTURE'))
print('Number of bands:', ds.RasterCount)

for i in range(1, ds.RasterCount+1):
    band = ds.GetRasterBand(i) # in GDAL, band are indexed starting at 1!
    interp = band.GetColorInterpretation()
    interp_name = gdal.GetColorInterpretationName(interp)
    (w,h)=band.GetBlockSize()
    print('Band %d, block size %dx%d, color interp %s' % (i,w,h,interp_name))

    ovr_count = band.GetOverviewCount()
    for j in range(ovr_count):
        ovr_band = band.GetOverview(j) # but overview bands starting at 0
        print(' Overview %d: %dx%d'%(j, ovr_band.XSize, ovr_band.YSize))

in_ds = gdal.Open(in_rst)
    in_gt = in_ds.GetGeoTransform()

    # TODO: Get number of bands (nBands) from input raster to pass on to driver.Create()

    inv_gt = gdal.InvGeoTransform(in_gt)
    if gdal.VersionInfo()[0] == '1':
        if inv_gt[0] == 1:
            inv_gt = inv_gt[1]
        else:
            raise RunTimeError('Inverse geotransform failed')
    elif inv_gt is None:
        raise RunTimeError('Inverse geotransform failed')

    # Compute upper left and lower right offsets
    offsets_ul = gdal.ApplyGeoTransform(inv_gt, aoi_ulx, aoi_uly)
    offsets_lr = gdal.ApplyGeoTransform(inv_gt, aoi_lrx, aoi_lry)
    off_ulx, off_uly = map(int, offsets_ul)
    off_lrx, off_lry = map(int, offsets_lr)

    # Compute number of rows and columns to extract
    rows = off_lry - off_uly
    columns = off_lrx - off_ulx

    gtiff_driver = gdal.GetDriverByName('GTiff')
    out_ds = gtiff_driver.Create(rst_out_fn, columns, rows, n_bands)
    out_ds.SetProjection(in_ds.GetProjection())

    # Put new origin coordinates in geotransform
    subset_ulx, subset_uly = gdal.ApplyGeoTransform(in_gt, off_ulx, off_uly)
    out_gt = list(in_gt)
    out_gt[0] = subset_ulx