Convert x,y coordinates into CRS:84 for GeoJSON

Question

I have several .tif (GeoTIFF) images each with different CRS and spatial resolutions. For each image (.tif file), I have the CRS, 4 extents, data type and pixel size. I then use these images in an annotation tool to annotate them. The annotation tool doesn't understand the geo-information from the .tif and the output annotations (XML files) from the tool have x,y coordinates.

How do I convert these x,y coordinates into CRS:84 coordinates required for the GeoJSON format? I have multiple such output XML files and would like to do it programmatically.

For example image 1 has the following properties,

Extent  367943.0000000000000000,6148320.0000000000000000 : 371030.0000000000000000,6150549.0000000000000000
Width   3087
Height  2229
Data type   Byte - Eight bit unsigned integer
Dimensions  X: 3087 Y: 2229 Bands: 3
Origin  367943,6.15055e+06
Pixel Size  1,-1

The annotated output of one polygon from the image has the following data,

<polygon label="Label name" occluded="0" source="manual" points="528.36,947.58;627.34,1101.98;730.14,1035.14;626.50,879.89" z_order="0">
      <attribute name="attriute name">Attribute value</attribute>
</polygon>

Each XML file has thousands of such polygons and I have hundreds of such XML files. The points here are x,y coordinates and I want to convert them into CRS:84 format.

Edit: In the original question i asked for help to convert it into EPSG:4326 but as @nmtoken pointed out GeoJSON uses CRS:84 (LONG/LAT) instead of EPSG:4326(LAT/LONG). I've updated the question accordingly.

Answer 1

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 coordinates\n{}'.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 coordinates\n{}'.format(a.ExportToWkt()))

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

    print('\nPolygon in CRS84 coordinates\n{}'.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.