[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