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

I have a source satellite image. It is a rectangular shape x by y. I also have the stored lon/lat coordinates of the corners and the center point in the xml meta file. Is there a way with the gdal or rasterio library in Python, to transform (scale, rotate, affine transform) the image to the shape defined in lon/lat coordinates of the corner points?

Answer

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(centerx)+"\n")
    worldfile.write(str(centery)+"\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]))

coo= 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]]