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

David M Nelson dnelson99 at gmail.com
Thu Jul 17 20:50:32 EDT 2008


Disregard the second question as we figured out that the second script
was generated by the gdalwarp of the first script.

Thanks again Frank!

Here's what we did:

Using 4 GCP for each of the image corners, we first derived the basic *.vrt.
The image filename and the vrt filename had to be the same for it to work.
All coordinates used were simple GCS.  The image was 16 bit single single
band.

<VRTDataset rasterXSize="16193" rasterYSize="16193">
 <Metadata/>
 <GCPList>
   <GCP Id="1" Pixel="0.0000" Line="0.0000" X="44.5962400" Y="17.7531250"/>
   <GCP Id="2" Pixel="16192.0000" Line="0.0000" X="50.4539309"
Y="15.4314191"/>
   <GCP Id="3" Pixel="16192.0000" Line="16192.0000" X="47.9521170"
Y="9.9314724"/>
   <GCP Id="4" Pixel="0.0000" Line="16192.0000" X="42.3048904"
Y="12.1465804"/>
 </GCPList>
 <VRTRasterBand dataType="UInt16" band="1">
   <Metadata/>
   <ColorInterp>Gray</ColorInterp>
   <SimpleSource>
     <SourceFilename relativeToVRT="1">LRG_IMAGE.tif</SourceFilename>
     <SourceBand>1</SourceBand>
     <SourceProperties RasterXSize="16193" RasterYSize="16193"
DataType="UInt16" BlockXSize="512" BlockYSize="512"/>
     <SrcRect xOff="0" yOff="0" xSize="16193" ySize="16193"/>
     <DstRect xOff="0" yOff="0" xSize="16193" ySize="16193"/>
   </SimpleSource>
 </VRTRasterBand>
</VRTDataset>

Then to create the output warped vrt, we simply performed the gdalwarp:
gdalwarp -of VRT LRG_IMAGE.vrt WARP2.vrt

The generated output file contained the following:

<VRTDataset rasterXSize="21926" rasterYSize="21045"
subClass="VRTWarpedDataset">
 <GeoTransform>  4.2252126314763458e+01,  3.7168429155331572e-04, 
0.0000000000000000e+00,  1.7726475525000001e+01,  0.0000000000000000e+00,
-3.7168429155331572e-04</GeoTransform>
  <VRTRasterBand dataType="UInt16" band="1" subClass="VRTWarpedRasterBand"/>
  <BlockXSize>512</BlockXSize>
  <BlockYSize>128</BlockYSize>
  <GDALWarpOptions>
    <WarpMemoryLimit>6.71089e+07</WarpMemoryLimit>
    <ResampleAlg>NearestNeighbour</ResampleAlg>
    <WorkingDataType>UInt16</WorkingDataType>
    <Option name="INIT_DEST">0</Option>
    <SourceDataset relativeToVRT="0">LRG_IMAGE.vrt</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="4.459624000000E+01" Y="1.775312500000E+01"/>
                  <GCP Id="2" Pixel="16192.0000" Line="0.0000"
X="5.045393090000E+01" Y="1.543141910000E+01"/>
                  <GCP Id="3" Pixel="16192.0000" Line="16192.0000"
X="4.795211700000E+01" Y="9.931472400000E+00"/>
                  <GCP Id="4" Pixel="0.0000" Line="16192.0000"
X="4.230489040000E+01" Y="1.214658040000E+01"/>
                </GCPList>
              </GCPTransformer>
            </SrcGCPTransformer>
           
<DstGeoTransform>42.25212631476346,0.0003716842915533157,0,17.726475525,0,-0.0003716842915533157</DstGeoTransform>
           
<DstInvGeoTransform>-113677.4603472922,2690.45537496587,0,47692.29135543721,0,-2690.45537496587</DstInvGeoTransform>
          </GenImgProjTransformer>
        </BaseTransformer>
      </ApproxTransformer>
    </Transformer>
    <BandList>
      <BandMapping src="1" dst="1"/>
    </BandList>
  </GDALWarpOptions>
</VRTDataset>

This script was readable in our application and projected the rotated image
onto the basemap.

Awesome!

David




David M Nelson wrote:
> 
> Thanks for the response, Frank.  I've been able to parse through and
> figure
> out most of the code.  I am a little unclear on the DstGeoTransform part
> of the second script.  Do I understand that the DstGeoTransform describes
> the original projection coordinates and the DstInvGeoTransform describes
> the destination projection coordinates (and that everything is in meters)?
> I am also not sure about the function of the GCP since the GeoTransform
> already defined the UL pixel value?  I just want to make sure that I
> understand the differences between the control point sets.
> 
> Thanks!
> 
> David
> 
> 
> 
> 
> 
> Frank Warmerdam-2 wrote:
>> 
>> 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
>> 
>> _______________________________________________
>> gdal-dev mailing list
>> gdal-dev at lists.osgeo.org
>> http://lists.osgeo.org/mailman/listinfo/gdal-dev
>> 
>> 
> 
> 

-- 
View this message in context: http://www.nabble.com/Using-VRT-to-project-rotate-tiff-files-on-the-fly-tp18448351p18520772.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.



More information about the gdal-dev mailing list