Python GDAL: Cómo obtener los errores dx/dy/residual para los GCPs tras el uso de gdalwarp

Question

Llevo unas semanas desarrollando mi plugin para QGIS. Va a ser bastante similar al conocido plugin GDAL Georeferencer, pero incluyendo algunas características más útiles. Hasta ahora mi código ha funcionado bastante bien, pero me he encontrado con un problema grave y no consigo resolverlo.

Me gustaría saber cómo obtener el dx/dy/residual ¿errores para los GCPs después de la excusión de gdalwarp?

Me refiero a los errores que muestra el georreferenciador GDAL en las 3 últimas columnas de la tabla GCP. Lo hicieron de alguna manera en C ++, pero no puedo averiguar cómo.

Lo que hago es simplemente:

1

gdal_translate -of GTiff  [GCPs] [source.tif] [temp.tif]

2

gdalwarp -r near -order 1 -co COMPRESS=NONE [temp.tif] [final.tif]

Entonces escribo gdalinfo [temp.tif] para que pueda ver la lista de GCP en la respuesta, pero es el archivo antes de la georeferencia. Desafortunadamente, gdalinfo [final.tif] no da ninguna información sobre los BPC.

Sin embargo, si utilizo VRT en lugar de GeoTIff para gdal_translate y gralwarp salidas hay un poco más de información. Creo que podría ser un paso más allá. En [final.vrt] Los GCPs están listados, pero sin errores (dx/dy/rms). También hay líneas denominadas <DstGeoTransform> y < DstInvGeoTransform >. Creo que pueden tener algo en común con los errores que quiero conseguir.

¿Alguien sabe cómo tratarla?

La solución perfecta para mí sería conseguir una fórmula para calcular las nuevas coordenadas de los GCPs utilizando la transformación de gdalwarp . Entonces simplemente compararía las coordenadas calculadas y las que el usuario recibió haciendo clic en el lienzo del mapa.

Desgraciadamente, gdalwarp trabaja con rasters, no con puntos puros.

¿Tiene alguna idea para resolver mi problema?

Ya he buscado en el stackexchange. El tema ya se ha planteado, pero no se ha resuelto. Cómo conservar los BPC después de gdalwarp

Adjunto una parte del ejemplo final.vrt debajo de

    <VRTDataset rasterXSize="547" rasterYSize="732" subClass="VRTWarpedDataset">
  <GeoTransform> 2.2187573162167482e+006, 3.0844837064488870e+000, 0.0000000000000000e+000, 6.4578332124060132e+006, 0.0000000000000000e+000,-3.0844837064488870e+000</GeoTransform>
  <VRTRasterBand dataType="Byte" band="1" subClass="VRTWarpedRasterBand">
    <ColorInterp>Palette</ColorInterp>
...
  </VRTRasterBand>
  <VRTRasterBand dataType="Byte" band="2" subClass="VRTWarpedRasterBand">
    <ColorInterp>Alpha</ColorInterp>
  </VRTRasterBand>
  <BlockXSize>512</BlockXSize>
  <BlockYSize>128</BlockYSize>
  <GDALWarpOptions>
    <WarpMemoryLimit>6.71089e+007</WarpMemoryLimit>
    <ResampleAlg>NearestNeighbour</ResampleAlg>
    <WorkingDataType>Byte</WorkingDataType>
    <Option name="INIT_DEST">0</Option>
    <SourceDataset relativeToVRT="0">G:Python/Python/okno_QGIS/smieci/mapa_krakowa 1.vrt</SourceDataset>
    <Transformer>
      <ApproxTransformer>
        <MaxError>0.125</MaxError>
        <BaseTransformer>
          <GenImgProjTransformer>
            <SrcGCPTransformer>
              <GCPTransformer>
                <Order>2</Order>
                <Reversed>0</Reversed>
                <GCPList>
                  <GCP Id="" Pixel="158.0590" Line="178.4810" X="2.219220000000E+006" Y="6.457240000000E+006" />
                  <GCP Id="" Pixel="374.7860" Line="297.1410" X="2.219910000000E+006" Y="6.456890000000E+006" />
                  <GCP Id="" Pixel="177.6730" Line="347.1550" X="2.219290000000E+006" Y="6.456710000000E+006" />
                  <GCP Id="" Pixel="79.6064" Line="301.0640" X="2.219010000000E+006" Y="6.456850000000E+006" />
                  <GCP Id="" Pixel="232.5900" Line="468.2670" X="2.219490000000E+006" Y="6.456350000000E+006" />
                  <GCP Id="" Pixel="63.4254" Line="426.0980" X="2.218950000000E+006" Y="6.456450000000E+006" />
                  <GCP Id="" Pixel="323.7910" Line="535.4420" X="2.219770000000E+006" Y="6.456150000000E+006" />
                  <GCP Id="" Pixel="191.8920" Line="616.8370" X="2.219360000000E+006" Y="6.455880000000E+006" />
                  <GCP Id="" Pixel="218.8600" Line="162.3000" X="2.219430000000E+006" Y="6.457290000000E+006" />
                  <GCP Id="" Pixel="303.1980" Line="201.5260" X="2.219690000000E+006" Y="6.457170000000E+006" />
                  <GCP Id="" Pixel="282.1130" Line="140.7250" X="2.219620000000E+006" Y="6.457360000000E+006" />
                  <GCP Id="" Pixel="107.0650" Line="102.4790" X="2.219080000000E+006" Y="6.457460000000E+006" />
                  <GCP Id="" Pixel="264.4610" Line="67.6657" X="2.219560000000E+006" Y="6.457600000000E+006" />
                  <GCP Id="" Pixel="400.2830" Line="117.1890" X="2.219990000000E+006" Y="6.457460000000E+006" />
                  <GCP Id="" Pixel="415.9740" Line="209.3720" X="2.220050000000E+006" Y="6.457160000000E+006" />
                </GCPList>
              </GCPTransformer>
            </SrcGCPTransformer>
            <DstGeoTransform>2218757.3162167482,3.084483706448887,0,6457833.2124060132,0,-3.084483706448887</DstGeoTransform>
            <DstInvGeoTransform>-719328.59025252087,0.32420336599906463,0,2093651.2645225821,0,-0.32420336599906463</DstInvGeoTransform>
          </GenImgProjTransformer>
        </BaseTransformer>
      </ApproxTransformer>
    </Transformer>
    <BandList>
      <BandMapping src="1" dst="1" />
    </BandList>
    <DstAlphaBand>2</DstAlphaBand>
  </GDALWarpOptions>
</VRTDataset>

Answer 1

Lamentablemente no he encontrado una solución fácil para mi problema, pero he escrito mi propio código que calcula lo que quiero.

Quizás a alguien le interese, así que he copiado aquí un poco de código. El código calcula los valores de error (dX, dY, dXY) para las transformaciones helicoidales/polinómicas basadas en los puntos de control del terreno (GCP). De hecho, el plugin GDAL Geo-referencer da resultados ligeramente diferentes en los mismos conjuntos de datos. Probablemente no utilizan el método de mínimos cuadrados que he elegido aquí, sino otra cosa. De todos modos, es un paso más allá, así que espero que mi código sea útil para alguien.

Si alguien sabe qué tipo de método (en lugar del método de mínimos cuadrados) utiliza GDAL Geo-referencer (el plugin de QGIS), estoy ansioso por saberlo.

import numpy as np
import math

#HELMERT  TRANSFORMATIONS
def helm_trans(gcps): #cgps is numpy.array [x, y, Xmap, Ymap]
    n = len(gcps)
    xo, yo, Xo, Yo = 0.0, 0.0, 0.0, 0.0

    #JANEK calculate center of gravity 
    for i in range(n): 
        xo, yo, Xo, Yo = xo + gcps[i,0], yo + gcps[i,1], Xo + gcps[i,2], Yo + gcps[i,3]
    xo, yo, Xo, Yo = xo/n, yo/n, Xo/n, Yo/n

    del_x, del_y, del_X, del_Y = gcps[:,0] - xo, gcps[:,1] - yo, gcps[:,2] - Xo, gcps[:,3] - Yo

    #JANEK calculation of unknowns
    a_up, a_down, b_up, b_down= 0, 0, 0, 0 
    for i in range(n):
        a_up += del_x[i]*del_Y[i] - del_y[i]*del_X[i]
        a_down += del_x[i]*del_x[i] + del_y[i]*del_y[i]
        b_up += del_x[i]*del_X[i] + del_y[i]*del_Y[i]

    b_down = a_down
    a = a_up/a_down
    b= b_up/b_down
    c = yo*a - xo*b + Xo
    d = -xo*a - yo*b + Yo # a,b,c,d are transformation parameters X = c + b*x - a*y / Y = d + a*x + b*y 

    #JANEK calculate new coordinates for points based on transformation values
    Xi = (gcps[:,0] - xo)*b - (gcps[:,1] - yo)*a + Xo 
    Yi = (gcps[:,0] - xo)*a + (gcps[:,1] - yo)*b + Yo

    #JANEK compare calculated values to "clicked" ones
    V_X = Xi - gcps[:,2] 
    V_Y = Yi - gcps[:,3]
    V_XY = np.sqrt(V_X*V_X + V_Y*V_Y)

    V_XY_sum_sq, V_X_sum_sq, V_Y_sum_sq = 0, 0, 0
    for i in range(n):
        V_XY_sum_sq += V_XY[i]*V_XY[i]
        V_X_sum_sq += V_X[i]*V_X[i]
        V_Y_sum_sq += V_Y[i]*V_Y[i]

    mo = math.sqrt(V_XY_sum_sq/(n)) #avarage error
    mox = math.sqrt(V_X_sum_sq/(n)) #avarage x error
    moy = math.sqrt(V_Y_sum_sq/(n)) #avarage y error

    return V_X, V_Y, V_XY, mo, mox, moy, [a, b, c, d]

#POLYNOMIAL TRANSFORMATIONS    
def polynomial(order, points_arr, Ax_row, Ay_row, LX_row, LY_row): #( {1, 2 or 3}, np.array, x-row, y-row, X-row, Y-row)
    n = len(points_arr)
    points = np.zeros((len(points_arr), 4), dtype=np.float)
    points[:, 0] = points_arr[:, Ax_row]
    points[:, 1] = points_arr[:, Ay_row]
    points[:, 2] = points_arr[:, LX_row]
    points[:, 3] = points_arr[:, LY_row]
    if order == 1:
        #X = a0 + a1x + a2y
        #Y = b0 + b1x + b2y
        Axy = np.zeros((len(points_arr), 3), dtype=np.float)
        Axy[:, 0] = 1
        Axy[:,1:3] = points[:, 0:2]
    elif order == 2:
        #X = a0 + a1x + a2y + a3xy + a4x^2 + a5y^2
        #Y = b0 + b1x + b2y + b3xy + b4x^2 + b5y^2
        Axy = np.zeros((len(points_arr), 6), dtype=np.float)
        Axy[:, 0] = 1 #a0
        Axy[:, 1] = points[ : , 0] # a1
        Axy[:, 2] = points[ : , 1] # a2
        Axy[:, 3] = points[ : , 0] * points[ : , 1] # ...
        Axy[:, 4] = points[ : , 0] * points[ : , 0]
        Axy[:, 5] = points[ : , 1] * points[ : , 1]
    elif order == 3:
        #X = a0 + a1x + a2y + a3xy + a4x^2 + a5y^2 + a6x^3 + a7x^2y + a8xy^2 + a9y^3
        #Y = b0 + b1x + b2y + b3xy + b4x^2 + b5y^2 + b6x^3 + b7x^2y + b8xy^2 + b9y^3
        Axy = np.zeros((len(points_arr), 10), dtype=np.float)
        Axy[:, 0] = 1 #a0
        Axy[:, 1] = points[ : , 0] # a1
        Axy[:, 2] = points[ : , 1] # a2
        Axy[:, 3] = points[ : , 0] * points[ : , 1] # ...
        Axy[:, 4] = points[ : , 0] * points[ : , 0]
        Axy[:, 5] = points[ : , 1] * points[ : , 1] #
        Axy[:, 6] = points[ : , 0] * points[ : , 0] * points[ : , 0]
        Axy[:, 7] = points[ : , 0] * points[ : , 0] * points[ : , 1]
        Axy[:, 8] = points[ : , 0] * points[ : , 1] * points[ : , 1]
        Axy[:, 9] = points[ : , 1] * points[ : , 1] * points[ : , 1]

    BX = points[ : , 2]
    BY = points[ : , 3]

    aaa_X = np.linalg.lstsq(Axy,BX) #transf parameters calculation using least square method
    bbb_Y = np.linalg.lstsq(Axy,BY)

    predXs = Axy.dot(aaa_X[0])
    predYs = Axy.dot(bbb_Y[0])

    V_X = points[ : , 2] - np.array(predXs)
    V_Y = points[ : , 3] - np.array(predYs)
    V_XY = np.sqrt(V_X*V_X + V_Y*V_Y)

    V_XY_sum_sq, V_X_sum_sq, V_Y_sum_sq = 0, 0, 0
    for i in range(n):
        V_XY_sum_sq += V_XY[i]*V_XY[i]
        V_X_sum_sq += V_X[i]*V_X[i]
        V_Y_sum_sq += V_Y[i]*V_Y[i]

    mo = math.sqrt(V_XY_sum_sq/(n)) #avarage error
    mox = math.sqrt(V_X_sum_sq/(n)) #avarage x error
    moy = math.sqrt(V_Y_sum_sq/(n)) #avarage y error

    return V_X, V_Y, V_XY, mo, mox, moy, aaa_X[0], bbb_Y[0]