1.. _rfc-22: 2 3================================================================================ 4RFC 22: RPC Georeferencing 5================================================================================ 6 7Author: Frank Warmerdam 8 9Contact: warmerdam@pobox.com 10 11Status: Adopted, Implemented 12 13Summary 14------- 15 16It is proposed that GDAL support an additional mechanism for geolocation 17of imagery based on rational polynomial coefficients (RPCs) represented 18as metadata. 19 20Many modern raw satellite products are distributed with RPCs, including 21products from GeoEye, and DigitalGlobe. RPCs provide a higher systematic 22description of georeferencing over an image, and also contain 23information on the viewing geometry that in theory makes orthocorrection 24(given a DEM) and some 3D operations like building height computation 25possible. 26 27RPC Domain Metadata 28------------------- 29 30Datasets with RPCs will include the following dataset level metadata 31items in the "RPC" domain to identify the rational polynomials. 32 33- ERR_BIAS: Error - Bias. The RMS bias error in meters per horizontal 34 axis of all points in the image (-1.0 if unknown) 35- ERR_RAND: Error - Random. RMS random error in meters per horizontal 36 axis of each point in the image (-1.0 if unknown) 37- LINE_OFF: Line Offset 38- SAMP_OFF: Sample Offset 39- LAT_OFF: Geodetic Latitude Offset 40- LONG_OFF: Geodetic Longitude Offset 41- HEIGHT_OFF: Geodetic Height Offset 42- LINE_SCALE: Line Scale 43- SAMP_SCALE: Sample Scale 44- LAT_SCALE: Geodetic Latitude Scale 45- LONG_SCALE: Geodetic Longitude Scale 46- HEIGHT_SCALE: Geodetic Height Scale 47- LINE_NUM_COEFF (1-20): Line Numerator Coefficients. Twenty 48 coefficients for the polynomial in the Numerator of the rn equation. 49 (space separated) 50- LINE_DEN_COEFF (1-20): Line Denominator Coefficients. Twenty 51 coefficients for the polynomial in the Denominator of the rn 52 equation. (space separated) 53- SAMP_NUM_COEFF (1-20): Sample Numerator Coefficients. Twenty 54 coefficients for the polynomial in the Numerator of the cn equation. 55 (space separated) 56- SAMP_DEN_COEFF (1-20): Sample Denominator Coefficients. Twenty 57 coefficients for the polynomial in the Denominator of the cn 58 equation. (space separated) 59 60These fields are directly derived from the document prospective GeoTIFF 61RPC document at: 62 63`http://geotiff.maptools.org/rpc_prop.html <http://geotiff.maptools.org/rpc_prop.html>`__ 64 65The line and pixel offset expressed with LINE_OFF and SAMP_OFF are with 66respect to the center of the pixel (#5993) 67 68Updating NITF Driver 69-------------------- 70 71- Already supports RPCs in this model, but will be modified to put them 72 in the RPC domain instead of the primary metadata domain. 73- Add support for reading Digital Globe .RPB files. 74- No support for writing RPCs for now. 75 76Updating GeoTIFF Driver 77----------------------- 78 79- Upgrade to support reading Digital Globe .RPB files. 80- Possible support reading Space Imaging (GeoEye?) rpc.txt files. 81- Support reading RPC TIFF tag (per 82 `http://geotiff.maptools.org/rpc_prop.html <http://geotiff.maptools.org/rpc_prop.html>`__) 83- Support writing RPC TIFF tag. 84- Support writing .RPB files (if RPB=YES or PROFILE not GDALGeoTIFF). 85 86Changes to GenImgProj Transformer 87--------------------------------- 88 89Currently it is difficult to reliably create a warp transformer based on 90RPCs using GDALGenImgProjTransformer() as it will use a geotransform in 91preference to RPCs if available. Many images with useful RPC information 92also include a geotransform (approximate or accurate). It is therefore 93proposed to modify the GDALCreateGenImgProjTransformer() function to 94make it practical to provide more direction in the creation of the 95transformer. The proposed new function is: 96 97:: 98 99 void * 100 GDALCreateGenImgProjTransformer2( GDALDatasetH hSrcDS, GDALDatasetH hDstDS, 101 char **papszOptions ); 102 103Supported Options: 104 105- SRC_SRS: WKT SRS to be used as an override for hSrcDS. 106- DST_SRS: WKT SRS to be used as an override for hDstDS. 107- GCPS_OK: If false, GCPs will not be used, default is TRUE. 108- MAX_GCP_ORDER: the maximum order to use for GCP derived polynomials 109 if possible. The default is to autoselect based on the number of 110 GCPs. A value of -1 triggers use of Thin Plate Spline instead of 111 polynomials. 112- METHOD: may have a value which is one of GEOTRANSFORM, 113 GCP_POLYNOMIAL, GCP_TPS, GEOLOC_ARRAY, RPC to force only one 114 geolocation method to be considered on the source dataset. 115- RPC_HEIGHT: A fixed height to be used with RPC calculations. 116 117This replaces the older function which did not include support for 118passing arbitrary options, and was thus not easily extended. The old 119function will be re-implemented with a call to the new functions. 120 121The most important addition is the METHOD option which can be set to 122specifically use one of the image to georeferenced coordinate system 123methods instead of leaving it up to the code to pick the one it thinks 124is best. 125 126Changes to gdalwarp and gdaltransform 127------------------------------------- 128 129In order to facilitate passing transformer options into the updated 130GDALCreateGenImgProjTransformer2(), the gdalwarp and gdaltransform 131programs (built on this function) will be updated to include a -to 132(transformer option) switch, and to use the new function. 133 134Preserving Geolocation Through Translation 135------------------------------------------ 136 137The RPC information needs to be copied and preserved through 138translations that do not alter the spatial arrangement of the data. To 139that end RPC metadata copying will be added to: 140 141- VRT driver's CreateCopy(). 142- GDALDriver's default CreateCopy(). 143- GDALPamDataset::CopyInfo() 144- gdal_translate will be updated to copy RPC metadata to the 145 intermediate internal VRT if, and only if, no resizing or subsetting 146 is being done. 147 148Changes to RPC Transformer 149-------------------------- 150 151- Implement iterative "back transform" from pixel/line to 152 lat/long/height instead of simple linear approximator. 153- Add support for RPC_HEIGHT offset, so all Z values to transformer are 154 assumed to be relative to this offset (normally really and average 155 elevation for the scene). 156- Make RPC Transformer serializable (in VRT files, etc). 157 158Backward Compatibility Issues 159----------------------------- 160 161Previously the NITF driver returned RPC metadata in the default domain. 162With the implementation of this RFC for GDAL 1.6.0 any applications 163using this metadata would need to consult the RPC domain instead. The 164RPC\_ prefix on the metadata values has also been removed. 165 166The GDALCreateGenImgProjTransformer() function is preserved, so no 167compatibility issues are anticipated by the addition of the new 168generalized factory. 169 170SWIG Bindings Issues 171-------------------- 172 173- The raw access is by the established metadata api, so no changes are 174 needed for this. 175- The Warp API is only bound at a high level, so there should be no 176 changes in this regard. 177- For testing purposes it is desirable to provide a binding around the 178 GDAL transformer API. The following planned binding is based loosely 179 on OGRCoordinateTransformation API binding. So far I have only found 180 the TransformPoint( bDstToSrc, x, y, z ) entry point to be useful in 181 Python and even that ends up returning a (bSuccess, (x, y, z)) result 182 which is somewhat awkward. Is there a better way of doing this? 183 184:: 185 186 /************************************************************************/ 187 /* Transformer */ 188 /************************************************************************/ 189 190 %rename (Transformer) GDALTransformerInfoShadow; 191 class GDALTransformerInfoShadow { 192 private: 193 GDALTransformerInfoShadow(); 194 public: 195 %extend { 196 197 GDALTransformerInfoShadow( GDALDatasetShadow *src, GDALDatasetShadow *dst, 198 char **options ) { 199 GDALTransformerInfoShadow *obj = (GDALTransformerInfoShadow*) 200 GDALCreateGenImgProjTransformer2( (GDALDatasetH)src, (GDALDatasetH)dst, 201 options ); 202 return obj; 203 } 204 205 ~GDALTransformerInfoShadow() { 206 GDALDestroyTransformer( self ); 207 } 208 209 // Need to apply argin typemap second so the numinputs=1 version gets applied 210 // instead of the numinputs=0 version from argout. 211 %apply (double argout[ANY]) {(double inout[3])}; 212 %apply (double argin[ANY]) {(double inout[3])}; 213 int TransformPoint( int bDstToSrc, double inout[3] ) { 214 int nRet, nSuccess = TRUE; 215 216 nRet = GDALUseTransformer( self, bDstToSrc, 217 1, &inout[0], &inout[1], &inout[2], 218 &nSuccess ); 219 220 return nRet && nSuccess; 221 } 222 %clear (double inout[3]); 223 224 int TransformPoint( double argout[3], int bDstToSrc, 225 double x, double y, double z = 0.0 ) { 226 int nRet, nSuccess = TRUE; 227 228 argout[0] = x; 229 argout[1] = y; 230 argout[2] = z; 231 nRet = GDALUseTransformer( self, bDstToSrc, 232 1, &argout[0], &argout[1], &argout[2], 233 &nSuccess ); 234 235 return nRet && nSuccess; 236 } 237 238 #ifdef SWIGCSHARP 239 %apply (double *inout) {(double*)}; 240 %apply (double *inout) {(int*)}; 241 #endif 242 int TransformPoints( int bDstToSrc, 243 int nCount, double *x, double *y, double *z, 244 int *panSuccess ) { 245 int nRet; 246 247 nRet = GDALUseTransformer( self, bDstToSrc, nCount, x, y, z, panSuccess ); 248 249 return nRet; 250 } 251 #ifdef SWIGCSHARP 252 %clear (double*); 253 %clear (int*); 254 #endif 255 256 } /*extend */ 257 }; 258 259Documentation 260------------- 261 262In addition to standard API documentation, the RPC metadata mechanism 263will be introduced into the "GDAL Data Model" document. 264 265Implementation 266-------------- 267 268This work will be implemented by Frank Warmerdam with support from the 269Canadian Nuclear Safety Commission. 270 271Testing 272------- 273 274- A test script for the transformer API covering RPC, GCP_TPS, 275 GCP_POLYNOMIAL, GEOLOC and GEOTRANSFORM methods will be implemented. 276- A test script for reading and writing RPB, and GeoTIFF RPC tags will 277 be written. 278