I'm trying to automate georeference of big amount of non/georeferenced image file stored in TIF format to work with them in QGIS 3.16.4-Hannover.
These TIF files are 8 bits (B/N) and I only have aux.xml sidecar file containing a list of GCP.
This is an example of a PAM Dataset *.aux.xml file:
<PAMDataset>
<SRS dataAxisToSRSAxisMapping="1,2">PROJCS["...",GEOGCS["...,"AUTHORITY["EPSG","23031"]]]</SRS>
<PAMRasterBand band="1">
<Metadata>
<MDI key="STATISTICS_APPROXIMATE">YES</MDI>
<MDI key="STATISTICS_MAXIMUM">250</MDI>
<MDI key="STATISTICS_MEAN">105.39136338541</MDI>
<MDI key="STATISTICS_MINIMUM">0</MDI>
<MDI key="STATISTICS_STDDEV">62.095857421551</MDI>
<MDI key="STATISTICS_VALID_PERCENT">100</MDI>
</Metadata>
</PAMRasterBand>
<GCPList>
<GCP Id="" Pixel="10627.14063" Line="2047.42188" X="303207.84625" Y="4709035.90307" />
<GCP Id="" Pixel="12700.51563" Line="13847.60938" X="303973.58269" Y="4705214.27235" />
<GCP Id="" Pixel="14836.50000" Line="9597.48438" X="304752.03808" Y="4706639.88411" />
... (more than 50 GCP points) ...
<GCP Id="" Pixel="9628.54688" Line="6340.40625" X="302924.73400" Y="4707741.75435" />
</GCPList>
</PAMDataset>
There are several possible transformation methods in GDAL [Ref1, pag 20]: TPS, Affine, 2nd and 3rd order polynomials. By default, polynomial order 1 will be selected if there are less than 6 points, otherwise polynomial order 2.
[Ref1] https://download.osgeo.org/gdal/workshop/foss4ge2015/workshop_gdal.pdf (Pg 20)
I have to maintain my list of GCP points but I need to use an affine transformation.
I have seen in some ESRI aux.xml file <PolynomialOrder>
tag, but QGIS doesn't recognize the <PolynomialOrder>1</PolynomialOrder>
instruction.
Is there any TAG in aux.xml file to set the transformation method to overpass GDAL default behaviour?
Best Answer
Use gdalwarp https://gdal.org/programs/gdalwarp.html with
-order 1
but instead of writing a new TIFF use VRT https://gdal.org/drivers/raster/vrt.html as outputformat.Command:
Open the image.vrt file with QGIS and it will apply the first order polynomial on-the-fly.
PolynomialOrder
appears only in the ESRI variant of aux.xml. GDAL does support also that but I do not know if it applies the PolynomialOrder instruction automatically. Maybe you can write ESRI style aux.xml by using https://github.com/OSGeo/gdal/blob/master/autotest/gcore/pam.py#L459 as a template and test.