[gdal-dev] vrt inverse RPC warping

[Even Rouault]

Hi Thomas,

(just seeing your message from the archive in 
https://lists.osgeo.org/pipermail/gdal-dev/2025-December/061216.html. 
There are issues with email delivering since last week with 
@lists.osgeo.org affecting a number of people)

I believe you can do what you want, but that will require some manual 
VRT editing (that could be automated in some script)

You need to first create a warped VRT from your raster, that you will 
use to copy/paste/update another warped VRT:

gdalwarp -r cubic -rpc raster.tif raster_rpc_warped.vrt

Then copy it into a dem_inverse_rpc_warped.vrt doing the following edits.

e.g in autotest we have

gdalwarp -r cubic -rpc autotest/alg/data/test_rpc_with_gt_bug_2460.tif 
in.vrt

the resulting in.vrt is the following. I've annoted the places that need 
to be modified to adapt it for your use case

<VRTDataset rasterXSize="281" rasterYSize="227" subClass="VRTWarpedDataset">
   <SRS dataAxisToSRSAxisMapping="2,1">GEOGCS["WGS 
84",DATUM["WGS_1984",SPHEROID["WGS 
84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]</SRS>

   <--- to be removed -->
   <GeoTransform> -1.1499575331273579e+02, 3.6616932415051353e-06,  
0.0000000000000000e+00, 3.6190393017469262e+01,  0.0000000000000000e+00, 
-3.6616932415051353e-06</GeoTransform>

   <VRTRasterBand dataType="Byte" band="1" 
subClass="VRTWarpedRasterBand">     <!--- update with your actual data 
type -->
     <ColorInterp>Gray</ColorInterp>
   </VRTRasterBand>
   <BlockXSize>281</BlockXSize>
   <BlockYSize>128</BlockYSize>
   <GDALWarpOptions>
     <WarpMemoryLimit>6.71089e+07</WarpMemoryLimit>
     <ResampleAlg>Cubic</ResampleAlg>
     <WorkingDataType>Byte</WorkingDataType>  <!--- update with your 
actual data type -->
     <Option name="INIT_DEST">0</Option>
     <Option name="SOURCE_EXTRA">5</Option>
     <Option name="ERROR_OUT_IF_EMPTY_SOURCE_WINDOW">FALSE</Option>
     <SourceDataset 
relativeToVRT="1">autotest/alg/data/test_rpc_with_gt_bug_2460.tif</SourceDataset> 
  <-- change here to your DEM
     <Transformer>
       <ApproxTransformer>
         <MaxError>0.125</MaxError>
         <BaseTransformer>
           <GenImgProjTransformer>
             <SrcRPCTransformer>                   <--- change this to 
be DstRPCTransformer
               <RPCTransformer>
                 <Reversed>0</Reversed>        <-- you may need to 
change that to 1. Not entirely sure ...
                 <HeightOffset>0</HeightOffset>
                 <PixErrThreshold>0.1</PixErrThreshold>
                 <Metadata> <--- to be replaced by what you get in in.vrt
                   <MDI key="ERR_BIAS">-1</MDI>
                   <MDI key="ERR_RAND">-1</MDI>
                   <MDI key="LINE_OFF">128</MDI>
                   <MDI key="SAMP_OFF">128</MDI>
                   <MDI key="LAT_OFF">36.1895</MDI>
                   <MDI key="LONG_OFF">-114.9943</MDI>
                   <MDI key="HEIGHT_OFF">500</MDI>
                   <MDI key="LINE_SCALE">128</MDI>
                   <MDI key="SAMP_SCALE">128</MDI>
                   <MDI key="LAT_SCALE">0.0015</MDI>
                   <MDI key="LONG_SCALE">0.0024</MDI>
                   <MDI key="HEIGHT_SCALE">1000</MDI>
                   <MDI key="MIN_LONG">-180</MDI>
                   <MDI key="MIN_LAT">-90</MDI>
                   <MDI key="MAX_LONG">180</MDI>
                   <MDI key="MAX_LAT">90</MDI>
                   <MDI key="LINE_NUM_COEFF">0.0418498012 -0.0149130525 
4.03738788 2.6824292 -0.293975583 -0.195230863 0.848327924 0.00111702848 
0.231318818 0.463079636 0.104397546 -7.92580025e-05 -0.0109232371 
0.073238365 0.0226941745 0.00203306277 0.138240173 0.0147348153 
-0.0944158497 0.133690662</MDI>
                   <MDI key="LINE_DEN_COEFF">0.966071333 -0.0705548972 
0.054937412 -0.143504966 -0.00263873607 0.0490624273 -0.0407232374 
0.00548420472 0.000493945772 -0.00508963274 0.00200894902 0.000300247695 
0.000485633376 -0.00767877615 -0.00041179685 -0.000543960426 
0.00619465736 -0.00264921176 -0.00147916846 -0.0158449621</MDI>
                   <MDI key="SAMP_NUM_COEFF">-0.0740525516 -5.21414237 
-0.00982810958 3.93746386 -0.262274148 -1.21406902 0.199674781 
0.416705428 -0.000130440488 0.673770884 0.104882676 -0.0562560276 
0.0158769672 -0.151393292 0.0156039821 3.67853087e-05 -0.0887577654 
-0.100124515 -0.0118096929 0.196201559</MDI>
                   <MDI key="SAMP_DEN_COEFF">0.965447163 -0.0775609504 
0.0482259778 -0.144208596 -0.0030041407 0.0508795964 -0.0371312046 
0.0105356006 -0.00287000921 -0.00501789141 0.00201318489 -7.63194545e-05 
0.00061804791 -0.00832315087 -0.000327509631 -0.000811519204 
0.0055600645 -0.00335719232 -0.000913845686 -0.0155829651</MDI>
                 </Metadata>
               </RPCTransformer>
             </SrcRPCTransformer>

             <-- you need to change this to SrcGeoTransform with the 
actual GeoTransform of your DEM -->
<DstGeoTransform>-114.99575331273579,3.6616932415051353e-06,0,36.190393017469262,0,-3.6616932415051353e-06</DstGeoTransform>

             <-- you need to change this to SrcInvGeoTransform with the 
actual inverse GeoTransform of your DEM -->

<DstInvGeoTransform>31405075.665340796,273097.69935532642,0,9883513.0718359239,0,-273097.69935532642</DstInvGeoTransform>
           </GenImgProjTransformer>
         </BaseTransformer>
       </ApproxTransformer>
     </Transformer>
     <BandList>
       <BandMapping src="1" dst="1" />
     </BandList>
   </GDALWarpOptions>
</VRTDataset>

Even

-- 
http://www.spatialys.com
My software is free, but my time generally not.