1

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.

4
  • 1
    Some additional information on how you create the XMLs would be helpful.
    – Erik
    Apr 12, 2021 at 12:57
  • 1
    into epsg:4326 coordinates required for the GeoJSON format? the CRS for GeoJSON is NOT epsg:4326, it's CRS:84. The difference is the order of the axes. EPSG:4326 is LAT/LONG whilst CRS:84 is LONG/LAT
    – nmtoken
    Apr 14, 2021 at 12:26
  • @Erik The XMLs are output automatically by the annotation software. There are a plethora of other data but this is the only part related to the coordinates. What other information would help? Apr 15, 2021 at 10:58
  • @nmtoken Thanks for this information. I looked into it and it does seem that you are correct. I'll update the original question accordingly. Apr 15, 2021 at 10:59

1 Answer 1

2

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.

5
  • Thanks for your response. I tried code and it seems to be working. I'm trying to understand what you're doing here, from what i understand you're getting the values of the origin from the original .tif file and scaling the x,y coordinates with pixel size and adding it to the origin. This should convert the x,y coordinates into the CRS of the input .tif file, but not necessarily into the required CRS:84 format correct? gdaltransform seems to be capable of performing this conversion but how do i programmatically set the -s_srs based on the input .tif? Apr 15, 2021 at 13:42
  • 1
    I am a little confused. You have some .tif rasters. Each one has an extent. You want to annotate this extent as lat-long, or do you want to annotate some other polygon inside the raster (as lat-long)?
    – wingnut
    Apr 15, 2021 at 14:22
  • I have some .tif rasters. I want to annotate this and then convert the output into GeoJSON format. I use some tool to help with this annotation, the output of this tool (annotated data) has x,y coordinate points and I want to convert them into CRS:84. As far as i understand your code helps convert this x,y coordinate from the output of the annotation back to the original CRS of the input image. But i want to convert this x,y to CRS:84 for further processing. Hopefully this clears up the confusion. Apr 15, 2021 at 14:32
  • I'll have a look tomorrow. If you are in a hurry, the basic approach is to transform your list of x,y for each annotation into a shapely Polygon, apply an affine transform to get the equivalent in raster CRS coordinates, then transform that to CRS84, and finally convert it to GeoJSON. You would then output that with the existing attributes in your XML annotations.
    – wingnut
    Apr 15, 2021 at 16:47
  • Modified the code so it is a function that takes the XML 'points' definition and generates a CRS84 JSON version.
    – wingnut
    Apr 16, 2021 at 8:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.