Convert x,y coordinates into CRS:84 for GeoJSON

The GDAL library (available in QGIS) can give you a GeoTransform from a file. Apply it to your x,y data. Example below.

Update

A function in a loop is probably what you are looking for.

import shapely.wkt as wkt
from shapely import affinity
from osgeo import gdal, ogr, osr

def transform_string(string,rasterfile):
    # make the transform parameters and affine transform matrix
    ds = gdal.Open(rasterfile)
    gt = ds.GetGeoTransform()
    gg = [gt[1],gt[2],gt[4],gt[5],gt[0],gt[3]]

    s_srs = ds.GetSpatialRef()
    # CRS84 definition
    aswkt = 'GEOGCRS["WGS 84 (CRS84)",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 longitude (Lon)",east,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic latitude (Lat)",north,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],USAGE[SCOPE["unknown"],AREA["World"],BBOX[-90,-180,90,180]],ID["OGC","CRS84"]]'
    t_srs = osr.SpatialReference(aswkt)
    
    # convert XML string to WKT polygon, still in pixel coordinates
    # The first point is repeated at the end to close the polygon.
    p1start = string.split(';')[0].replace(',' ,' ')
    p1 = 'POLYGON (('+string.replace(',',' ').replace(';',',')+','+p1start+'))'
    raster_poly = wkt.loads(p1)
    
    print('Polygon in pixel coordinatesn{}'.format(raster_poly.wkt))
    
    # apply affine transform to make the data compatible with s_srs 
    a = ogr.CreateGeometryFromWkt(affinity.affine_transform(raster_poly,gg).wkt,s_srs)
    
    print('nPolygon in local coordinatesn{}'.format(a.ExportToWkt()))

    # transform from s_srs to CRS84
    a.TransformTo(t_srs)

    print('nPolygon in CRS84 coordinatesn{}'.format(a.ExportToWkt()))

    return a.ExportToJson()


transform_string("528.36,947.58;627.34,1101.98;730.14,1035.14;626.50,879.89",'l4.tif')

Running this gives me:

Polygon in raster coordinates
POLYGON ((528.36 947.58, 627.34 1101.98, 730.14 1035.14, 626.5 879.89, 528.36 947.58))

Polygon in local coordinates
POLYGON ((611696.73980244 1884952.55916432,650781.02476286 1791845.54610192,691373.71514406 1832151.71678256,650449.3335185 1925771.30083656,611696.73980244 1884952.55916432))

Polygon in CRS84 coordinates
POLYGON ((150.444464266743 -73.1017810420974,151.879033180361 -73.9098595628977,153.052208429374 -73.5169246485067,151.540097611741 -72.7136192802267,150.444464266743 -73.1017810420974))

The returned JSON looks like this:

'{ "type": "Polygon", "coordinates": [ [ [ 150.444464266742955, -73.101781042097386 ], [ 151.879033180361176, -73.909859562897665 ], [ 153.052208429373934, -73.516924648506688 ], [ 151.540097611741402, -72.713619280226709 ], [ 150.444464266742955, -73.101781042097386 ] ] ] }'

Hopefully, this can be modified to give you what you need.

Note:

1. The basic process is from pixel to local to OGC:CRS84 coordinates

2. The input CRS has been determined from the GeoTiff, to avoid any errors in working off a separate file. It is reordered to transform the polygon in GDAL, which seems to use a different order of coefficients.

3. The CRS of the output layer (t_srs) is also determined in the function, based on CRS84's WKT definition. This should be moved outside the function and passed in as a parameter as it only needs to be done once.

4. CRS 84 will probably look the same as EPSG:4326

5. You could pass in a list of text-formatted polygon definitions, and loop through them, returning a list of JSON (or other format) outputs.

6. Each XML polygon needs to have the first point repeated in order to form a closed polygon.

7. I used an arbitrary GeoTiff file for testing, so the coordinates won't be the same as yours.

Tested Code in Python Console

from osgeo import gdal
import os
hp = QgsProject.instance().homePath()
fname = os.path.join(hp,'l4.tif')
print(fname)
ds = gdal.Open(fname)
print(ds)
gt = ds.GetGeoTransform()
print(gt)
x = 10
y = 10
a = gdal.ApplyGeoTransform(gt, x,y)
print(a)

You can get the inverse transform from:

inv_gt = gdal.InvGeoTransform(gt)

if you need it.