1

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?

1
  • I know qgis does this when opening the tif file. Jun 7 '19 at 14:15
3

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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