Transformación de una imagen a una forma definida por las esquinas de la imagen en coordenadas terrestres - Python

Question

Tengo una imagen de satélite de origen. Es una forma rectangular x por y. También tengo las coordenadas lon/lat almacenadas de las esquinas y el punto central en el archivo meta xml. ¿Hay alguna manera con el gdal o rasterio en Python, para transformar (escalar, rotar, transformar afinadamente) la imagen a la forma definida en coordenadas lon/lat de los puntos de esquina?

Answer 1

En otras palabras, quiere crear un Archivo mundial a partir de las coordenadas de las 4 esquinas y de la anchura y altura de la imagen

1. Se obtiene la anchura y la altura de la imagen con osgeo.gdal , rasterio o cualquier otra biblioteca para abrir archivos de imagen como Almohada y otros.

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

2. necesitas extraer los valores x e y de las esquinas del archivo XML (no conozco la estructura del archivo XML)

   Mis esquinas:

   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. crear un Archivo mundial con un simple script (sin gdal ni rasterio) con sólo 2 puntos de esquina

   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. con gdal y los 4 puntos de esquina como puntos de control del terreno.

   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. Hay otras soluciones como tab2tfw.py de Michael Kalbermatten (es exactamente el mismo problema que el archivo de fichas de MapInfo) o utilizando affine6p , nudged-py o Ajuste_afín para estimar los parámetros de transformación afín entre dos conjuntos de puntos 2D, pero hay que tener en cuenta que los datos ráster, procedentes de sus orígenes de procesamiento de imágenes, utilizan un sistema de referencia diferente para acceder a los píxeles (véase Transformación afín en Python )

Ejemplo con Numpy y los 4 puntos de esquina (transformación Affine) ( 0,0 origen en la parte superior izquierda )

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

Ejemplo de resultado con nugged ( 0,0 origen en la parte superior izquierda )

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