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

[Hays Barrett]

That works great! Thanks Evan! 


From: "Even Rouault" <even.rouault at spatialys.com> 
To: gdal-dev at lists.osgeo.org 
Cc: "Hays Barrett" <hbarrett at edac.unm.edu> 
Sent: Tuesday, April 21, 2020 7:18:04 AM 
Subject: Re: [gdal-dev] Help with finding residual error for a given GCP. 



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/20200422/94f51f88/attachment-0001.html>