1

I have hundreds of very large TIFFs which were georeferenced 15+ years ago using ArcGIS 9.2. Viewing them in ArcMap 9.2 the images are properly located and sized/warped on the map. However, they need to be rectified to the final projection. The GCPs are embedded in the TIF and also in the tif.xml file.
While I can achieve warping to a final rectified image using the Clip tool, it takes a long time for each image. I've tried using GDAL warp utilities to do this, but am getting empty final images (all black).
eg:

\OSGeo4W\bin\gdalwarp.exe -overwrite -co COMPRESS=DEFLATE IQALUIT_1948_A11535_043.tif out.tif

\OSGeo4W\bin\gdalwarp.exe -overwrite --config GDAL_GEOREF_SOURCES PAM -co COMPRESS=DEFLATE IQALUIT_1948_A11535_043.tif out.tif

Here is the output from \OSGeo4W\bin\gdalinfo.exe IQALUIT_1948_A11535_043.tif:

Driver: GTiff/GeoTIFF
Files: IQALUIT_1948_A11535_043.tif
       IQALUIT_1948_A11535_043.aux
       IQALUIT_1948_A11535_043.rrd
       IQALUIT_1948_A11535_043.tif.aux.xml
Size is 18565, 18492
GCP Projection =
PROJCRS["Canada_Lambert_Conformal_Conic",
    BASEGEOGCRS["NAD83",
        DATUM["North American Datum 1983",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6269]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Lambert Conic Conformal (2SP)",
            ID["EPSG",9802]],
        PARAMETER["Latitude of false origin",40,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8821]],
        PARAMETER["Longitude of false origin",-96,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8822]],
        PARAMETER["Latitude of 1st standard parallel",50,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Latitude of 2nd standard parallel",70,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8824]],
        PARAMETER["Easting at false origin",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8826]],
        PARAMETER["Northing at false origin",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8827]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]
Data axis to CRS axis mapping: 1,2
GCP[  0]: Id=, Info=
          (16.5886813769484,11.2421187314794) -> (1298609.77217853,2919147.39696529,0)
GCP[  1]: Id=, Info=
          (16.8939942179049,11.7402607351451) -> (1298500.0219128,2919153.17329506,0)
GCP[  2]: Id=, Info=
          (18.8921073940182,10.0052797183835) -> (1298631.40807544,2919617.77480046,0)
GCP[  3]: Id=, Info=
          (21.4572641379127,13.0572806764721) -> (1297892.90562009,2919783.3167529,0)
GCP[  4]: Id=, Info=
          (19.6528779424261,13.8720869249719) -> (1297910.23024105,2919417.7417536,0)
GCP[  5]: Id=, Info=
          (20.6880569996538,12.0845177405408) -> (1298127.44460173,2919739.01694754,0)
GCP[  6]: Id=, Info=
          (21.5855745124354,11.8220205578855) -> (1298091.07382506,2919911.77813673,0)
GCP[  7]: Id=, Info=
          (19.6699549399509,11.5899398962065) -> (1298297.65636365,2919614.6948765,0)
GCP[  8]: Id=, Info=
          (19.2596578670837,14.6946775183915) -> (1297815.43048519,2919292.29625073,0)
GCP[  9]: Id=, Info=
          (18.8790359460127,12.7118694177316) -> (1298182.37712172,2919377.79059004,0)
GCP[ 10]: Id=, Info=
          (17.0335864905755,13.3239667836351) -> (1298226.2344595,2919045.4523473,0)
Metadata:
  TIFFTAG_MAXSAMPLEVALUE=255
  TIFFTAG_MINSAMPLEVALUE=0
  TIFFTAG_RESOLUTIONUNIT=3 (pixels/cm)
  TIFFTAG_XRESOLUTION=800.26886
  TIFFTAG_YRESOLUTION=800.14081
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,18492.0)
Upper Right (18565.0,    0.0)
Lower Right (18565.0,18492.0)
Center      ( 9282.5, 9246.0)
Band 1 Block=18565x1 Type=Byte, ColorInterp=Gray
  Overviews: 4642x4623, 2321x2312, 1161x1156, 581x578, 291x289, 146x145, 73x73, 37x37

and the IQALUIT_1948_A11535_043.tif.aux.xml file:

<?xml version="1.0" encoding="utf-8" ?>
<GeodataXform xsi:type='typens:PolynomialXform' xmlns:xsi='http://www.w3.org/2001/XMLSchema-instance'
              xmlns:xs='http://www.w3.org/2001/XMLSchema' xmlns:typens='http://www.esri.com/schemas/ArcGIS/9.2'>
    <PolynomialOrder>1</PolynomialOrder>
    <SpatialReference xsi:type='typens:ProjectedCoordinateSystem'>
        <WKT>PROJCS[&quot;Canada_Lambert_Conformal_Conic&quot;,GEOGCS[&quot;GCS_North_American_1983&quot;,DATUM[&quot;D_North_American_1983&quot;,SPHEROID[&quot;GRS_1980&quot;,6378137.0,298.257222101]],PRIMEM[&quot;Greenwich&quot;,0.0],UNIT[&quot;Degree&quot;,0.0174532925199433]],PROJECTION[&quot;Lambert_Conformal_Conic&quot;],PARAMETER[&quot;False_Easting&quot;,0.0],PARAMETER[&quot;False_Northing&quot;,0.0],PARAMETER[&quot;Central_Meridian&quot;,-96.0],PARAMETER[&quot;Standard_Parallel_1&quot;,50.0],PARAMETER[&quot;Standard_Parallel_2&quot;,70.0],PARAMETER[&quot;Latitude_Of_Origin&quot;,40.0],UNIT[&quot;Meter&quot;,1.0]]</WKT>
        <XOrigin>-35105400</XOrigin>
        <YOrigin>-29255500</YOrigin>
        <XYScale>128287945.084531</XYScale>
        <ZOrigin>-100000</ZOrigin>
        <ZScale>10000</ZScale>
        <MOrigin>-100000</MOrigin>
        <MScale>10000</MScale>
        <XYTolerance>1.55899293474704E-08</XYTolerance>
        <ZTolerance>0.0002</ZTolerance>
        <MTolerance>0.0002</MTolerance>
        <HighPrecision>true</HighPrecision>
    </SpatialReference>
    <SourceGCPs xsi:type='typens:ArrayOfDouble'>
        <Double>16.5886813769484</Double>
        <Double>11.2421187314794</Double>
        <Double>16.8939942179049</Double>
        <Double>11.7402607351451</Double>
        <Double>18.8921073940182</Double>
        <Double>10.0052797183835</Double>
        <Double>21.4572641379127</Double>
        <Double>13.0572806764721</Double>
        <Double>19.6528779424261</Double>
        <Double>13.8720869249719</Double>
        <Double>20.6880569996538</Double>
        <Double>12.0845177405408</Double>
        <Double>21.5855745124354</Double>
        <Double>11.8220205578855</Double>
        <Double>19.6699549399509</Double>
        <Double>11.5899398962065</Double>
        <Double>19.2596578670837</Double>
        <Double>14.6946775183915</Double>
        <Double>18.8790359460127</Double>
        <Double>12.7118694177316</Double>
        <Double>17.0335864905755</Double>
        <Double>13.3239667836351</Double>
    </SourceGCPs>
    <TargetGCPs xsi:type='typens:ArrayOfDouble'>
        <Double>1298609.77217853</Double>
        <Double>2919147.39696529</Double>
        <Double>1298500.0219128</Double>
        <Double>2919153.17329506</Double>
        <Double>1298631.40807544</Double>
        <Double>2919617.77480046</Double>
        <Double>1297892.90562009</Double>
        <Double>2919783.3167529</Double>
        <Double>1297910.23024105</Double>
        <Double>2919417.7417536</Double>
        <Double>1298127.44460173</Double>
        <Double>2919739.01694754</Double>
        <Double>1298091.07382506</Double>
        <Double>2919911.77813673</Double>
        <Double>1298297.65636365</Double>
        <Double>2919614.6948765</Double>
        <Double>1297815.43048519</Double>
        <Double>2919292.29625073</Double>
        <Double>1298182.37712172</Double>
        <Double>2919377.79059004</Double>
        <Double>1298226.2344595</Double>
        <Double>2919045.4523473</Double>
    </TargetGCPs>
</GeodataXform>

I don't understand the units of the ESRI sourceGCPs - why the values are so small. Is it an offset instead of a pixel position, or coordinate?

What is the correct gdalwarp command to use to get this to rectify properly? (the same projected CRS as the source can be used for the target.)

(I've already looked at Create rectified GeoTIFF given aux.xml produced from ArcGIS)

Improve this question
3
  • How does the gdalinfo about out.tif look like?
    – user30184
    Commented Mar 20, 2023 at 7:07
  • Can you share the data? GDAL should be able to convert the small source pixel/line values from the aux.xml file into full pixels but with your dataset the conversion fails for some reason.
    – user30184
    Commented Mar 20, 2023 at 9:05
  • Follow the progress of github.com/OSGeo/gdal/issues/7484. ESRI has a special alternative for giving the source coordinates, if I understand right it means centimeters or inches from the corner if the original was a paper print.
    – user30184
    Commented Mar 21, 2023 at 7:56

1 Answer 1

Reset to default
2

GDAL did not support just that combination of ESRI variant of ground control points, TIFF, and aux.xml file. GDAL developers have verified and fixed the issue https://github.com/OSGeo/gdal/issues/7484. Install the development version of GDAL or build it from the sources and you can warp your images correctly.

Improve this answer
1
  • Thanks. I'll give it a go.
    – Matt
    Commented Mar 22, 2023 at 18:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.