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