[Gdal-dev] Using VRT to project/rotate tiff files on-the-fly

[Frank Warmerdam]

David M Nelson wrote:
> I have a series of unprojected, single-frame images.  We want to be able to
> project them on-the-fly by generating a VRT file for each image.  However,
> none of the images are oriented north-up--they are all uniquely oriented
> because of the flightpath taken by the platform during the time of
> acquisition.  I have not been able to find a way to generate a VRT that
> takes orientation into account.  I tried to use the X-rotation and the
> Y-rotation, but this completely distorts the image setting some of the
> coordinate values to 35,000+ deg.  I think that the orientation values are
> supposed to be for the reference frame anyway and not the images (?).
> 
> So my question is: is there a way to generate a VRT file which will project
> and orient correctly an unprojected image?

David,

If they are still very regularly sampled (only requiring an affine
transformation to a north up uniform grid) you should be able to express
the georeferencing with an appropriate geotransform though computing it may
be a bit tricky.

If the are not necessarily really regular you could create a VRT that
expresses their georeferencing in terms of georeferenced locations at
known pixel control points (ie. GCPs). A simple example of an image
georeferenced with GCPs using a VRT is:

<VRTDataset rasterXSize="512" rasterYSize="512">
   <Metadata/>
   <GCPList>
     <GCP Id="1" Pixel="0.0000" Line="0.0000" X="1.000000000000E+04" 
Y="5.000000000000E+03"/>
     <GCP Id="2" Pixel="500.0000" Line="0.0000" X="1.200000000000E+04" 
Y="7.000000000000E+03"/>
     <GCP Id="3" Pixel="250.0000" Line="500.0000" X="9.000000000000E+03" 
Y="6.000000000000E+03"/>
   </GCPList>
   <VRTRasterBand dataType="Byte" band="1">
     <Metadata/>
     <ColorInterp>Gray</ColorInterp>
     <SimpleSource>
       <SourceFilename relativeToVRT="1">gcp.tif</SourceFilename>
       <SourceBand>1</SourceBand>
       <SourceProperties RasterXSize="512" RasterYSize="512" DataType="Byte" 
BlockXSize="512" BlockYSize="16"/>
       <SrcRect xOff="0" yOff="0" xSize="512" ySize="512"/>
       <DstRect xOff="0" yOff="0" xSize="512" ySize="512"/>
     </SimpleSource>
   </VRTRasterBand>
</VRTDataset>

However, this depends on the application reading the image knowing how
to utilize GCPs which many do not.  You could actually rectify the image
with gdalwarp.  But if you don't want to do that you can also create a
"warped vrt" which is a virtual image representing what gdalwarp would
have produced.  For instance, for the above this might look like:

<VRTDataset rasterXSize="1448" rasterYSize="724" subClass="VRTWarpedDataset">
   <GeoTransform>  7.9520000000000000e+03,  2.8284271247461898e+00, 
0.0000000000000000e+00,  7.0480000000000000e+03,  0.0000000000000000e+00, 
-2.8284271247461898e+00</GeoTransform>
   <VRTRasterBand dataType="Byte" band="1" subClass="VRTWarpedRasterBand"/>
   <BlockXSize>512</BlockXSize>
   <BlockYSize>128</BlockYSize>
   <GDALWarpOptions>
     <WarpMemoryLimit>6.71089e+07</WarpMemoryLimit>
     <ResampleAlg>NearestNeighbour</ResampleAlg>
     <WorkingDataType>Byte</WorkingDataType>
     <Option name="INIT_DEST">0</Option>
     <SourceDataset relativeToVRT="1">gcp.tif</SourceDataset>
     <Transformer>
       <ApproxTransformer>
         <MaxError>0.125</MaxError>
         <BaseTransformer>
           <GenImgProjTransformer>
             <SrcGCPTransformer>
               <GCPTransformer>
                 <Order>1</Order>
                 <Reversed>0</Reversed>
                 <GCPList>
                   <GCP Id="1" Pixel="0.0000" Line="0.0000" 
X="1.000000000000E+04" Y="5.000000000000E+03"/>
                   <GCP Id="2" Pixel="500.0000" Line="0.0000" 
X="1.200000000000E+04" Y="7.000000000000E+03"/>
                   <GCP Id="3" Pixel="250.0000" Line="500.0000" 
X="9.000000000000E+03" Y="6.000000000000E+03"/>
                 </GCPList>
               </GCPTransformer>
             </SrcGCPTransformer>
 
<DstGeoTransform>7952,2.82842712474619,0,7048,0,-2.82842712474619</DstGeoTransform>
 
<DstInvGeoTransform>-2811.456561997713,0.3535533905932738,0,2491.844296901394,0,-0.3535533905932738</DstInvGeoTransform>
           </GenImgProjTransformer>
         </BaseTransformer>
       </ApproxTransformer>
     </Transformer>
     <BandList>
       <BandMapping src="1" dst="1"/>
     </BandList>
   </GDALWarpOptions>
</VRTDataset>

Any GDAL supporting application will read the warped vrt as a north up image.
The commands I used to create the above examples are:

gdal_translate plain.tif gcp.tif -gcp 0 0 10000 5000
         -gcp 500 0 12000 7000 -gcp 250 500 9000 6000

gdal_translate gcp.tif gcp.vrt -of VRT

gdalwarp -of VRT gcp.tif warped.vrt

Best regards,
-- 
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | President OSGeo, http://osgeo.org