Python GDAL: How to obtain dx/dy/residual errors for GCPs after gdalwarp usage

Question

I've been developing my QGIS plugin for few weeks now. It's going to be pretty similar to well known GDAL Georeferencer plugin, but including few more helpful features. So far my code has been working pretty fine, but I met a serious problem and I can't manage to solve it.

I'd love to know how to obtain the dx/dy/residual errors for GCPs after gdalwarp excuting?

I mean the errors that GDAL georeferencer shows in last 3 columns of GCP table. They did it somehow in C++, but I can't figure out how.

What I do is simply:

1

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

2

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

Then I type gdalinfo [temp.tif] so I can see GCP list in answer, but it is the file before georeference. Unfortunatlly, gdalinfo [final.tif] gives is no info about GCPs.

However if I use VRT format instead of GeoTIff for gdal_translate and gralwarp outputs there is a little bit more information. I think it could be a step beyond. In [final.vrt] GCPs are listed, but without errors (dx/dy/rms). There are also lines called <DstGeoTransform> and <DstInvGeoTransform>. I think they may have something in common with errors I want to get.

Does anybody know how to deal with it?

The perfect solution for me would be to get a formula to calculate new coordinates for GCPs using the transformation from gdalwarp. Then i'd simply compare computed coordinates and those which user received clicking in map canvas.

Unfortunately, gdalwarp works with rasters, not with pure points.

Have you got any ideas to solve my problem?

I've already searched through the stackexchange. The subect have already been raised, but not solved. How to preserve GCP's after gdalwarp

I attach a part of example final.vrt below

    <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

Unfortunately I haven't found an easy solution for my problem, but I've written my own code that computes what I want.

Maybe someone will be interested in it, so I copied a bit of code here. The code computes error values (dX, dY, dXY) for helmert/polynomial transformations based on Ground Control Points (GCP). In fact, the GDAL Geo-referencer plugin gives slightly different results on the same data-sets. Probably they don't use the least squares method which I've chosen here, but something else. Anyway, it is a step beyond, so hopefully my code will be helpful to somebody.

If somebody knows what kind of method (instead of least square method) is used by GDAL Geo-referencer (the QGIS plugin), I'm anxious to hear it.

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]