1 /******************************************************************************
2 *
3 * Project: Mapinfo Image Warper
4 * Purpose: Commandline program for doing a variety of image warps, including
5 * image reprojection.
6 * Author: Frank Warmerdam <warmerdam@pobox.com>
7 *
8 ******************************************************************************
9 * Copyright (c) 2002, i3 - information integration and imaging
10 * Fort Collin, CO
11 *
12 * Permission is hereby granted, free of charge, to any person obtaining a
13 * copy of this software and associated documentation files (the "Software"),
14 * to deal in the Software without restriction, including without limitation
15 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 * and/or sell copies of the Software, and to permit persons to whom the
17 * Software is furnished to do so, subject to the following conditions:
18 *
19 * The above copyright notice and this permission notice shall be included
20 * in all copies or substantial portions of the Software.
21 *
22 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28 * DEALINGS IN THE SOFTWARE.
29 ****************************************************************************/
30
31 #include "gdal_alg.h"
32 #include "cpl_string.h"
33 #include "ogr_srs_api.h"
34
35 CPL_CVSID("$Id: gdalwarpsimple.c ff8146d84de7cba8e09d212d5481ea7d2ede3e98 2017-06-27 20:47:31Z Even Rouault $")
36
37 static GDALDatasetH
38 GDALWarpCreateOutput( GDALDatasetH hSrcDS, const char *pszFilename,
39 const char *pszFormat, const char *pszSourceSRS,
40 const char *pszTargetSRS, int nOrder,
41 char **papszCreateOptions );
42
43 static double dfMinX=0.0, dfMinY=0.0, dfMaxX=0.0, dfMaxY=0.0;
44 static double dfXRes=0.0, dfYRes=0.0;
45 static int nForcePixels=0, nForceLines=0;
46
47 /************************************************************************/
48 /* Usage() */
49 /************************************************************************/
50
Usage()51 static void Usage()
52
53 {
54 printf(
55 "Usage: gdalwarpsimple [--version] [--formats]\n"
56 " [-s_srs srs_def] [-t_srs srs_def] [-order n] [-et err_threshold]\n"
57 " [-te xmin ymin xmax ymax] [-tr xres yres] [-ts width height]\n"
58 " [-of format] [-co \"NAME=VALUE\"]* srcfile dstfile\n" );
59 exit( 1 );
60 }
61
62 /************************************************************************/
63 /* SanitizeSRS */
64 /************************************************************************/
65
SanitizeSRS(const char * pszUserInput)66 char *SanitizeSRS( const char *pszUserInput )
67
68 {
69 OGRSpatialReferenceH hSRS;
70 char *pszResult = NULL;
71
72 CPLErrorReset();
73
74 hSRS = OSRNewSpatialReference( NULL );
75 if( OSRSetFromUserInput( hSRS, pszUserInput ) == OGRERR_NONE )
76 OSRExportToWkt( hSRS, &pszResult );
77 else
78 {
79 CPLError( CE_Failure, CPLE_AppDefined,
80 "Translating source or target SRS failed:\n%s",
81 pszUserInput );
82 exit( 1 );
83 }
84
85 OSRDestroySpatialReference( hSRS );
86
87 return pszResult;
88 }
89
90 /************************************************************************/
91 /* main() */
92 /************************************************************************/
93
main(int argc,char ** argv)94 int main( int argc, char ** argv )
95
96 {
97 GDALDatasetH hSrcDS, hDstDS;
98 const char *pszFormat = "GTiff";
99 char *pszTargetSRS = NULL;
100 char *pszSourceSRS = NULL;
101 const char *pszSrcFilename = NULL, *pszDstFilename = NULL;
102 int bCreateOutput = FALSE, i, nOrder = 0;
103 void *hTransformArg, *hGenImgProjArg=NULL, *hApproxArg=NULL;
104 char **papszWarpOptions = NULL;
105 double dfErrorThreshold = 0.125;
106 GDALTransformerFunc pfnTransformer = NULL;
107 char **papszCreateOptions = NULL;
108
109 GDALAllRegister();
110
111 /* -------------------------------------------------------------------- */
112 /* Parse arguments. */
113 /* -------------------------------------------------------------------- */
114 for( i = 1; i < argc; i++ )
115 {
116 if( EQUAL(argv[i],"--version") )
117 {
118 printf( "%s\n", GDALVersionInfo( "--version" ) );
119 exit( 0 );
120 }
121 else if( EQUAL(argv[i],"--formats") )
122 {
123 int iDr;
124
125 printf( "Supported Formats:\n" );
126 for( iDr = 0; iDr < GDALGetDriverCount(); iDr++ )
127 {
128 GDALDriverH hDriver = GDALGetDriver(iDr);
129
130 printf( " %s: %s\n",
131 GDALGetDriverShortName( hDriver ),
132 GDALGetDriverLongName( hDriver ) );
133 }
134
135 exit( 0 );
136 }
137 else if( EQUAL(argv[i],"-co") && i < argc-1 )
138 {
139 papszCreateOptions = CSLAddString( papszCreateOptions, argv[++i] );
140 bCreateOutput = TRUE;
141 }
142 else if( EQUAL(argv[i],"-of") && i < argc-1 )
143 {
144 pszFormat = argv[++i];
145 bCreateOutput = TRUE;
146 }
147 else if( EQUAL(argv[i],"-t_srs") && i < argc-1 )
148 {
149 pszTargetSRS = SanitizeSRS(argv[++i]);
150 }
151 else if( EQUAL(argv[i],"-s_srs") && i < argc-1 )
152 {
153 pszSourceSRS = SanitizeSRS(argv[++i]);
154 }
155 else if( EQUAL(argv[i],"-order") && i < argc-1 )
156 {
157 nOrder = atoi(argv[++i]);
158 }
159 else if( EQUAL(argv[i],"-et") && i < argc-1 )
160 {
161 dfErrorThreshold = CPLAtof(argv[++i]);
162 }
163 else if( EQUAL(argv[i],"-tr") && i < argc-2 )
164 {
165 dfXRes = CPLAtof(argv[++i]);
166 dfYRes = fabs(CPLAtof(argv[++i]));
167 bCreateOutput = TRUE;
168 }
169 else if( EQUAL(argv[i],"-ts") && i < argc-2 )
170 {
171 nForcePixels = atoi(argv[++i]);
172 nForceLines = atoi(argv[++i]);
173 bCreateOutput = TRUE;
174 }
175 else if( EQUAL(argv[i],"-te") && i < argc-4 )
176 {
177 dfMinX = CPLAtof(argv[++i]);
178 dfMinY = CPLAtof(argv[++i]);
179 dfMaxX = CPLAtof(argv[++i]);
180 dfMaxY = CPLAtof(argv[++i]);
181 bCreateOutput = TRUE;
182 }
183 else if( argv[i][0] == '-' )
184 Usage();
185 else if( pszSrcFilename == NULL )
186 pszSrcFilename = argv[i];
187 else if( pszDstFilename == NULL )
188 pszDstFilename = argv[i];
189 else
190 Usage();
191 }
192
193 if( pszDstFilename == NULL )
194 Usage();
195
196 /* -------------------------------------------------------------------- */
197 /* Open source dataset. */
198 /* -------------------------------------------------------------------- */
199 hSrcDS = GDALOpen( pszSrcFilename, GA_ReadOnly );
200
201 if( hSrcDS == NULL )
202 exit( 2 );
203
204 /* -------------------------------------------------------------------- */
205 /* Check that there's at least one raster band */
206 /* -------------------------------------------------------------------- */
207 if ( GDALGetRasterCount(hSrcDS) == 0 )
208 {
209 fprintf(stderr, "Input file %s has no raster bands.\n", pszSrcFilename );
210 exit( 2 );
211 }
212
213 if( pszSourceSRS == NULL )
214 {
215 if( GDALGetProjectionRef( hSrcDS ) != NULL
216 && strlen(GDALGetProjectionRef( hSrcDS )) > 0 )
217 pszSourceSRS = CPLStrdup(GDALGetProjectionRef( hSrcDS ));
218
219 else if( GDALGetGCPProjection( hSrcDS ) != NULL
220 && strlen(GDALGetGCPProjection(hSrcDS)) > 0
221 && GDALGetGCPCount( hSrcDS ) > 1 )
222 pszSourceSRS = CPLStrdup(GDALGetGCPProjection( hSrcDS ));
223 else
224 pszSourceSRS = CPLStrdup("");
225 }
226
227 if( pszTargetSRS == NULL )
228 pszTargetSRS = CPLStrdup(pszSourceSRS);
229
230 /* -------------------------------------------------------------------- */
231 /* Does the output dataset already exist? */
232 /* -------------------------------------------------------------------- */
233 CPLPushErrorHandler( CPLQuietErrorHandler );
234 hDstDS = GDALOpen( pszDstFilename, GA_Update );
235 CPLPopErrorHandler();
236
237 if( hDstDS != NULL && bCreateOutput )
238 {
239 fprintf( stderr,
240 "Output dataset %s exists,\n"
241 "but some commandline options were provided indicating a new dataset\n"
242 "should be created. Please delete existing dataset and run again.",
243 pszDstFilename );
244 exit( 1 );
245 }
246
247 /* -------------------------------------------------------------------- */
248 /* If not, we need to create it. */
249 /* -------------------------------------------------------------------- */
250 if( hDstDS == NULL )
251 {
252 hDstDS = GDALWarpCreateOutput( hSrcDS, pszDstFilename, pszFormat,
253 pszSourceSRS, pszTargetSRS, nOrder,
254 papszCreateOptions );
255 papszWarpOptions = CSLSetNameValue( papszWarpOptions, "INIT", "0" );
256 CSLDestroy( papszCreateOptions );
257 papszCreateOptions = NULL;
258 }
259
260 if( hDstDS == NULL )
261 exit( 1 );
262
263 /* -------------------------------------------------------------------- */
264 /* Create a transformation object from the source to */
265 /* destination coordinate system. */
266 /* -------------------------------------------------------------------- */
267 hTransformArg = hGenImgProjArg =
268 GDALCreateGenImgProjTransformer( hSrcDS, pszSourceSRS,
269 hDstDS, pszTargetSRS,
270 TRUE, 1000.0, nOrder );
271
272 if( hTransformArg == NULL )
273 exit( 1 );
274
275 pfnTransformer = GDALGenImgProjTransform;
276
277 /* -------------------------------------------------------------------- */
278 /* Warp the transformer with a linear approximator unless the */
279 /* acceptable error is zero. */
280 /* -------------------------------------------------------------------- */
281 if( dfErrorThreshold != 0.0 )
282 {
283 hTransformArg = hApproxArg =
284 GDALCreateApproxTransformer( GDALGenImgProjTransform,
285 hGenImgProjArg, dfErrorThreshold );
286 pfnTransformer = GDALApproxTransform;
287 }
288
289 /* -------------------------------------------------------------------- */
290 /* Now actually invoke the warper to do the work. */
291 /* -------------------------------------------------------------------- */
292 GDALSimpleImageWarp( hSrcDS, hDstDS, 0, NULL,
293 pfnTransformer, hTransformArg,
294 GDALTermProgress, NULL, papszWarpOptions );
295
296 CSLDestroy( papszWarpOptions );
297
298 if( hApproxArg != NULL )
299 GDALDestroyApproxTransformer( hApproxArg );
300
301 if( hGenImgProjArg != NULL )
302 GDALDestroyGenImgProjTransformer( hGenImgProjArg );
303
304 /* -------------------------------------------------------------------- */
305 /* Cleanup. */
306 /* -------------------------------------------------------------------- */
307 GDALClose( hDstDS );
308 GDALClose( hSrcDS );
309
310 GDALDumpOpenDatasets( stderr );
311
312 GDALDestroyDriverManager();
313
314 exit( 0 );
315 }
316
317 /************************************************************************/
318 /* GDALWarpCreateOutput() */
319 /* */
320 /* Create the output file based on various commandline options, */
321 /* and the input file. */
322 /************************************************************************/
323
324 static GDALDatasetH
GDALWarpCreateOutput(GDALDatasetH hSrcDS,const char * pszFilename,const char * pszFormat,const char * pszSourceSRS,const char * pszTargetSRS,int nOrder,char ** papszCreateOptions)325 GDALWarpCreateOutput( GDALDatasetH hSrcDS, const char *pszFilename,
326 const char *pszFormat, const char *pszSourceSRS,
327 const char *pszTargetSRS, int nOrder,
328 char **papszCreateOptions )
329
330 {
331 GDALDriverH hDriver;
332 GDALDatasetH hDstDS;
333 void *hTransformArg;
334 double adfDstGeoTransform[6];
335 int nPixels=0, nLines=0;
336 GDALColorTableH hCT;
337
338 /* -------------------------------------------------------------------- */
339 /* Find the output driver. */
340 /* -------------------------------------------------------------------- */
341 hDriver = GDALGetDriverByName( pszFormat );
342 if( hDriver == NULL
343 || GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL ) == NULL )
344 {
345 int iDr;
346
347 printf( "Output driver `%s' not recognised or does not support\n",
348 pszFormat );
349 printf( "direct output file creation. The following format drivers are configured\n"
350 "and support direct output:\n" );
351
352 for( iDr = 0; iDr < GDALGetDriverCount(); iDr++ )
353 {
354 GDALDriverH hDriver = GDALGetDriver(iDr);
355
356 if( GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL) != NULL )
357 {
358 printf( " %s: %s\n",
359 GDALGetDriverShortName( hDriver ),
360 GDALGetDriverLongName( hDriver ) );
361 }
362 }
363 printf( "\n" );
364 exit( 1 );
365 }
366
367 /* -------------------------------------------------------------------- */
368 /* Create a transformation object from the source to */
369 /* destination coordinate system. */
370 /* -------------------------------------------------------------------- */
371 hTransformArg =
372 GDALCreateGenImgProjTransformer( hSrcDS, pszSourceSRS,
373 NULL, pszTargetSRS,
374 TRUE, 1000.0, nOrder );
375
376 if( hTransformArg == NULL )
377 return NULL;
378
379 /* -------------------------------------------------------------------- */
380 /* Get approximate output definition. */
381 /* -------------------------------------------------------------------- */
382 if( GDALSuggestedWarpOutput( hSrcDS,
383 GDALGenImgProjTransform, hTransformArg,
384 adfDstGeoTransform, &nPixels, &nLines )
385 != CE_None )
386 return NULL;
387
388 GDALDestroyGenImgProjTransformer( hTransformArg );
389
390 /* -------------------------------------------------------------------- */
391 /* Did the user override some parameters? */
392 /* -------------------------------------------------------------------- */
393 if( dfXRes != 0.0 && dfYRes != 0.0 )
394 {
395 CPLAssert( nPixels == 0 && nLines == 0 );
396 if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 )
397 {
398 dfMinX = adfDstGeoTransform[0];
399 dfMaxX = adfDstGeoTransform[0] + adfDstGeoTransform[1] * nPixels;
400 dfMaxY = adfDstGeoTransform[3];
401 dfMinY = adfDstGeoTransform[3] + adfDstGeoTransform[5] * nLines;
402 }
403
404 nPixels = (int) ((dfMaxX - dfMinX + (dfXRes/2.0)) / dfXRes);
405 nLines = (int) ((dfMaxY - dfMinY + (dfYRes/2.0)) / dfYRes);
406 adfDstGeoTransform[0] = dfMinX;
407 adfDstGeoTransform[3] = dfMaxY;
408 adfDstGeoTransform[1] = dfXRes;
409 adfDstGeoTransform[5] = -dfYRes;
410 }
411
412 else if( nForcePixels != 0 && nForceLines != 0 )
413 {
414 if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 )
415 {
416 dfMinX = adfDstGeoTransform[0];
417 dfMaxX = adfDstGeoTransform[0] + adfDstGeoTransform[1] * nPixels;
418 dfMaxY = adfDstGeoTransform[3];
419 dfMinY = adfDstGeoTransform[3] + adfDstGeoTransform[5] * nLines;
420 }
421
422 dfXRes = (dfMaxX - dfMinX) / nForcePixels;
423 dfYRes = (dfMaxY - dfMinY) / nForceLines;
424
425 adfDstGeoTransform[0] = dfMinX;
426 adfDstGeoTransform[3] = dfMaxY;
427 adfDstGeoTransform[1] = dfXRes;
428 adfDstGeoTransform[5] = -dfYRes;
429
430 nPixels = nForcePixels;
431 nLines = nForceLines;
432 }
433
434 else if( dfMinX != 0.0 || dfMinY != 0.0 || dfMaxX != 0.0 || dfMaxY != 0.0 )
435 {
436 dfXRes = adfDstGeoTransform[1];
437 dfYRes = fabs(adfDstGeoTransform[5]);
438
439 nPixels = (int) ((dfMaxX - dfMinX + (dfXRes/2.0)) / dfXRes);
440 nLines = (int) ((dfMaxY - dfMinY + (dfYRes/2.0)) / dfYRes);
441
442 adfDstGeoTransform[0] = dfMinX;
443 adfDstGeoTransform[3] = dfMaxY;
444 }
445
446 /* -------------------------------------------------------------------- */
447 /* Create the output file. */
448 /* -------------------------------------------------------------------- */
449 printf( "Creating output file is that %dP x %dL.\n", nPixels, nLines );
450
451 hDstDS = GDALCreate( hDriver, pszFilename, nPixels, nLines,
452 GDALGetRasterCount(hSrcDS),
453 GDALGetRasterDataType(GDALGetRasterBand(hSrcDS,1)),
454 papszCreateOptions );
455
456 if( hDstDS == NULL )
457 return NULL;
458
459 /* -------------------------------------------------------------------- */
460 /* Write out the projection definition. */
461 /* -------------------------------------------------------------------- */
462 GDALSetProjection( hDstDS, pszTargetSRS );
463 GDALSetGeoTransform( hDstDS, adfDstGeoTransform );
464
465 /* -------------------------------------------------------------------- */
466 /* Copy the color table, if required. */
467 /* -------------------------------------------------------------------- */
468 hCT = GDALGetRasterColorTable( GDALGetRasterBand(hSrcDS,1) );
469 if( hCT != NULL )
470 GDALSetRasterColorTable( GDALGetRasterBand(hDstDS,1), hCT );
471
472 return hDstDS;
473 }
474