[gdal-dev] Help with finding residual error for a given GCP.

[Even Rouault]

Hays,

GCPsToGeoTransform() is not used directly by gdalwarp (although it should likely be 
equivalent to running -order 1 polynomial), but by other code (OZI .map exporting, .tab 
exporting)

You rather want to use the gdal.Transformer() API which is tested in
https://github.com/OSGeo/gdal/blob/master/autotest/gcore/transformer.py

Even


> Hello all,
> I want to get the error for each GCP before warping so I can calculate the
> RMSE before georeferencing. With a list of GCPs, is it possible to get the
> eventual map location of each pixel coordinate to calculate the
> transformation error for each GCP? I have tried to use GCPsToGeoTransform
> but that gave me a map location that was further from the intended chosen
> GCP coordinate than the actual post- georeferencing pixel location.
> 
> Here is my GCPsToGeoTransform code:
> ##CODE_START
> 
> from osgeo import gdal
> #list of GCPs
> gcpjson={'0': {'row': 5305.626973704321, 'col': 1158.432032601927, 'lat':
> -107.18221364936802, 'lon': 35.08107313264506}, '1': {'row':
> 778.3230540174817, 'col': 2317.0737231100225, 'lat': -107.17634116783263,
> 'lon': 35.09958594602314}, '2': {'row': 3950.8777704175645, 'col':
> 4880.465723118029, 'lat': -107.163460244414, 'lon': 35.08679518563491},
> '3': {'row': 7609.882307742781, 'col': 4407.027704042081, 'lat':
> -107.16570574683065, 'lon': 35.07164256912238}, '4': {'row':
> 1822.0904204049455, 'col': 6486.900377637105, 'lat': -107.15529960157271,
> 'lon': 35.096105471961536}, '5': {'row': 6291.396402113635, 'col':
> 7216.1033738731185, 'lat': -107.15180240538085, 'lon': 35.077317337370516}}
> 
> #New list to put in gdal.GCP objects into
> gcpList=[]
> #Load GCP objects into gcpList
> for gcp in gcpjson:
> gcpList.append(gdal.GCP(gcpjson[gcp]['lat'],gcpjson[gcp]['lon'],0,gcpjson[gc
> p]['col'],gcpjson[gcp]['row']))
> 
> #create GeoTranform from gcp list
> GT = gdal.GCPsToGeoTransform(tuple(gcpList))
> 
> #Get pixel location from first GCP as a test.
> col=gcpjson['0']['col']
> row=gcpjson['0']['row']
> xp = GT[0] + col*GT[1] + row*GT[2]
> yp = GT[3] + col*GT[4] + row*GT[5]
> print("POINT ("+str(xp)+" "+str(yp)+" 0)")
> 
> ##CODE_END
> 
> This prints out the WKT 'POINT (-107.18216298 35.0808978834 0)' which is
> kinda close to -107.18221364936802, 35.08107313264506 but looking at my
> post warped image the chosen pixel is right on top of the intended
> coordinate so I am clearly doing something (or a lot) wrong.
> 
> 
> 
> Any help would be great.
> Thanks,
> -Hays


-- 
Spatialys - Geospatial professional services
http://www.spatialys.com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20200421/c9239f94/attachment-0001.html>