Transforming a image to a shape defined by image corners in earth coordinates - Python

In other words, you want to create a World file from the coordinates of the 4 corners and the width and height of the image

1. You get the width and height of the image with osgeo.gdal, rasterio or any other libraries to open image files as Pillow and others.

   dataset = rasterio.open('satel.tif')
   rasterx = dataset.width
   rastery = dataset.height
   

2. you need to extract the x and y values of the corners from the XML file (I don't know the structure of the XML file)

   My corners:

   upper left  (162013.302, 138172.271) 
   lower left  (162013.302, 128171.074) 
   upper right (170015.863, 138172.271) 
   lower right (170015.863, 128171.074) 
   

3. create a World file with a simple script (without gdal or rasterio) with only 2 corners points

   x1,y1 = [162013.302,  138172.271] # upper left corner
   x2,y2 = [170015.863,  128171.074  # lower right corner
   rasterx  = 7988
   rastery =  9983
   # pixel size 
   #x-component of the pixel width (x scale)
   ppx = (x2-x1)/rasterx
   # y-component of the pixel height (y scale)
   ppy = (y2-y1)/rastery
   print(ppx, ppy)
   (1.0018228592889353, -1.0018227987578898)
   # x-coordinate of the center of the original image's upper left pixel transformed to the map
   xcenter = x1 + (ppx * .5)
   # y-coordinate of the center of the original image's upper left pixel transformed to the map
   ycenter = y1 + (ppy * .5)
   print(xcenter,ycenter)
   (162013.80291142964, 138171.77008860064)
   # write the worldfile
   with open('satel.tfw', "w") as worldfile:
       worldfile.write(str(ppx)+"n")
       worldfile.write(str(0)+"n") # y-component of the pixel width (y-skew), generaly = 0 
       worldfile.write(str(0)+"n") # y-component of the pixel width (y-skew), generaly = 0 
       worldfile.write(str(ppy)+"n")
       worldfile.write(str(xcenter)+"n")
       worldfile.write(str(ycenter)+"n")
   

4. with gdal and the 4 corners points as ground control points.

   from osgeo import gdal
   fp= [[0,rasterx,rasterx,0],[0,0,rastery,rastery]]
   tp= [[162013.302, 170015.863, 170015.863, 162013.302], [ 128171.074,128171.074,   138172.271,   138172.271]]
   pix = list(zip(fp[0],fp[1]))
   coor= list(zip(tp[0],tp[1]))
   # compute the gdal.GCP parameters
   gcps = []
   for index, txt in enumerate(pix):
       gcps.append(gdal.GCP())
       gcps[index].GCPPixel = pix[index][0]
       gcps[index].GCPLine = rastery-int(pix[index][1])
       gcps[index].GCPX = coor[index][0]
       gcps[index].GCPY = coor[index][1]
       geotransform = gdal.GCPsToGeoTransform( gcps )
       print(geotransform)
       (162013.302, 1.0018228592889353, 0.0, 138172.271, 0.0, -1.0018227987578898)
       xcenter = geotransform[0] + (geotransform[1] * .5)
       ycenter = geotransform[3] + (geotransform[5] * .5)
       print(xcenter,ycenter)
       (162013.80291142964, 138171.77008860064)
       # write the worldfile
       ...
   

5. There are other solutions like tab2tfw.py of Michael Kalbermatten (this is exactly the same problem as MapInfo tab file) or using affine6p, nudged-py or Affine_Fit to estimate the affine transformation parameters between two sets of 2D points but but be careful that raster data, coming from its image processing origins, uses a different referencing system to access pixels (see Python affine transforms)

Example with Numpy and the 4 corners points (Affine transformation) ( 0,0 origin in the upper left )

import numpy as np
fp = np.matrix([[0,rasterx,rasterx,0],[0,0,rastery,rastery]])
newline = [1,1,1,1]
fp  = np.vstack([fp,newline])
tp = np.matrix([[162013.302, 170015.863, 170015.863, 162013.302], [ 128171.074,128171.074,   138172.271,   138172.271]])
M = tp * fp.I
print(M)
matrix([[     1.00182286,      0.        , 162013.302     ],
        [     0.        ,      1.0018228 , 128171.074     ]])

Example of result with nugged ( 0,0 origin in the upper left )

import nudged
from_pt = [(0, 0), (rasterx, 0), (rasterx, rastery), (0, rastery)]
to_pt = [(162013.302, 128171.074), (170015.863, 128171.074), (170015.863, 138172.271), (162013.302, 138172.271)]
trans = nudged.estimate(from_pt, to_pt)
print(trans.get_matrix())
[[1.0018228223855314,0, 162013.3021473922],        
[0, 1.0018228223855314,   128171.0738820626], 
[0, 0, 1]]