source: trunk/zoo-project/zoo-services/gdal/warp/service.c

Last change on this file was 369, checked in by djay, 7 years ago

Add GCP option to gdal_translate service. Add warp service based on gdal_warp.

File size: 81.0 KB
Line 
1/******************************************************************************
2 * $Id: gdalwarp.cpp 24214 2012-04-08 20:17:17Z etourigny $
3 *
4 * Project:  High Performance Image Reprojector
5 * Purpose:  Test program for high performance warper API.
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 "gdalwarper.h"
32#include "cpl_string.h"
33#include "ogr_spatialref.h"
34#include "ogr_api.h"
35#include <vector>
36#include "service.h"
37
38CPL_CVSID("$Id: gdalwarp.cpp 24214 2012-04-08 20:17:17Z etourigny $");
39
40extern "C"{
41static void
42LoadCutline( const char *pszCutlineDSName, const char *pszCLayer, 
43             const char *pszCWHERE, const char *pszCSQL, 
44             void **phCutlineRet );
45static void
46TransformCutlineToSource( GDALDatasetH hSrcDS, void *hCutline,
47                          char ***ppapszWarpOptions, char **papszTO );
48
49static GDALDatasetH
50GDALWarpCreateOutput( maps*& conf,char **papszSrcFiles, const char *pszFilename, 
51                      const char *pszFormat, char **papszTO,
52                      char ***ppapszCreateOptions, GDALDataType eDT,
53                      void ** phTransformArg,
54                      GDALDatasetH* phSrcDS );
55
56static double          dfMinX=0.0, dfMinY=0.0, dfMaxX=0.0, dfMaxY=0.0;
57static double          dfXRes=0.0, dfYRes=0.0;
58static int             bTargetAlignedPixels = FALSE;
59static int             nForcePixels=0, nForceLines=0, bQuiet = FALSE;
60static int             bEnableDstAlpha = FALSE, bEnableSrcAlpha = FALSE;
61
62static int             bVRT = FALSE;
63
64/******************************************************************************/
65/*! \page gdalwarp gdalwarp
66
67image reprojection and warping utility
68
69\section gdalwarp_synopsis SYNOPSIS
70
71\htmlonly
72Usage:
73\endhtmlonly
74
75\verbatim
76gdalwarp [--help-general] [--formats]
77    [-s_srs srs_def] [-t_srs srs_def] [-to "NAME=VALUE"]
78    [-order n | -tps | -rpc | -geoloc] [-et err_threshold]
79    [-refine_gcps tolerance [minimum_gcps]]
80    [-te xmin ymin xmax ymax] [-tr xres yres] [-tap] [-ts width height]
81    [-wo "NAME=VALUE"] [-ot Byte/Int16/...] [-wt Byte/Int16]
82    [-srcnodata "value [value...]"] [-dstnodata "value [value...]"] -dstalpha
83    [-r resampling_method] [-wm memory_in_mb] [-multi] [-q]
84    [-cutline datasource] [-cl layer] [-cwhere expression]
85    [-csql statement] [-cblend dist_in_pixels] [-crop_to_cutline]
86    [-of format] [-co "NAME=VALUE"]* [-overwrite]
87    srcfile* dstfile
88\endverbatim
89
90\section gdalwarp_description DESCRIPTION
91
92<p>
93The gdalwarp utility is an image mosaicing, reprojection and warping
94utility. The program can reproject to any supported projection,
95and can also apply GCPs stored with the image if the image is "raw"
96with control information.
97
98<p>
99<dl>
100<dt> <b>-s_srs</b> <em>srs def</em>:</dt><dd> source spatial reference set.
101The coordinate systems that can be passed are anything supported by the
102OGRSpatialReference.SetFromUserInput() call, which includes EPSG PCS and GCSes
103(ie. EPSG:4296), PROJ.4 declarations (as above), or the name of a .prf file
104containing well known text.</dd>
105<dt> <b>-t_srs</b> <em>srs_def</em>:</dt><dd> target spatial reference set.
106The coordinate systems that can be passed are anything supported by the
107OGRSpatialReference.SetFromUserInput() call, which includes EPSG PCS and GCSes
108(ie. EPSG:4296), PROJ.4 declarations (as above), or the name of a .prf file
109containing well known text.</dd>
110<dt> <b>-to</b> <em>NAME=VALUE</em>:</dt><dd> set a transformer option suitable
111to pass to GDALCreateGenImgProjTransformer2(). </dd>
112<dt> <b>-order</b> <em>n</em>:</dt><dd> order of polynomial used for warping
113(1 to 3). The default is to select a polynomial order based on the number of
114GCPs.</dd>
115<dt> <b>-tps</b>:</dt><dd>Force use of thin plate spline transformer based on
116available GCPs.</dd>
117<dt> <b>-rpc</b>:</dt> <dd>Force use of RPCs.</dd>
118<dt> <b>-geoloc</b>:</dt><dd>Force use of Geolocation Arrays.</dd>
119<dt> <b>-et</b> <em>err_threshold</em>:</dt><dd> error threshold for
120transformation approximation (in pixel units - defaults to 0.125).</dd>
121<dt> <b>-refine_gcps</b> <em>tolerance minimum_gcps</em>:</dt><dd>  (GDAL >= 1.9.0) refines the GCPs by automatically eliminating outliers.
122Outliers will be eliminated until minimum_gcps are left or when no outliers can be detected.
123The tolerance is passed to adjust when a GCP will be eliminated.
124Not that GCP refinement only works with polynomial interpolation.
125The tolerance is in pixel units if no projection is available, otherwise it is in SRS units.
126If minimum_gcps is not provided, the minimum GCPs according to the polynomial model is used.</dd>
127<dt> <b>-te</b> <em>xmin ymin xmax ymax</em>:</dt><dd> set georeferenced
128extents of output file to be created (in target SRS).</dd>
129<dt> <b>-tr</b> <em>xres yres</em>:</dt><dd> set output file resolution (in
130target georeferenced units)</dd>
131<dt> <b>-tap</b>:</dt><dd> (GDAL >= 1.8.0) (target aligned pixels) align
132the coordinates of the extent of the output file to the values of the -tr,
133such that the aligned extent includes the minimum extent.</dd>
134<dt> <b>-ts</b> <em>width height</em>:</dt><dd> set output file size in
135pixels and lines. If width or height is set to 0, the other dimension will be
136guessed from the computed resolution. Note that -ts cannot be used with -tr</dd>
137<dt> <b>-wo</b> <em>"NAME=VALUE"</em>:</dt><dd> Set a warp options.  The
138GDALWarpOptions::papszWarpOptions docs show all options.  Multiple
139 <b>-wo</b> options may be listed.</dd>
140<dt> <b>-ot</b> <em>type</em>:</dt><dd> For the output bands to be of the
141indicated data type.</dd>
142<dt> <b>-wt</b> <em>type</em>:</dt><dd> Working pixel data type. The data type
143of pixels in the source image and destination image buffers.</dd>
144<dt> <b>-srcnodata</b> <em>value [value...]</em>:</dt><dd> Set nodata masking
145values for input bands (different values can be supplied for each band).  If
146more than one value is supplied all values should be quoted to keep them
147together as a single operating system argument.  Masked values will not be
148used in interpolation.  Use a value of <tt>None</tt> to ignore intrinsic nodata settings on the source dataset.</dd>
149<dt> <b>-dstnodata</b> <em>value [value...]</em>:</dt><dd> Set nodata values
150for output bands (different values can be supplied for each band).  If more
151than one value is supplied all values should be quoted to keep them together
152as a single operating system argument.  New files will be initialized to this
153value and if possible the nodata value will be recorded in the output
154file.</dd>
155<dt> <b>-dstalpha</b>:</dt><dd> Create an output alpha band to identify
156nodata (unset/transparent) pixels. </dd>
157<dt> <b>-wm</b> <em>memory_in_mb</em>:</dt><dd> Set the amount of memory (in
158megabytes) that the warp API is allowed to use for caching.</dd>
159<dt> <b>-multi</b>:</dt><dd> Use multithreaded warping implementation.
160Multiple threads will be used to process chunks of image and perform
161input/output operation simultaneously.</dd>
162<dt> <b>-q</b>:</dt><dd> Be quiet.</dd>
163<dt> <b>-of</b> <em>format</em>:</dt><dd> Select the output format. The default is GeoTIFF (GTiff). Use the short format name. </dd>
164<dt> <b>-co</b> <em>"NAME=VALUE"</em>:</dt><dd> passes a creation option to
165the output format driver. Multiple <b>-co</b> options may be listed. See
166format specific documentation for legal creation options for each format.
167</dd>
168
169<dt> <b>-cutline</b> <em>datasource</em>:</dt><dd>Enable use of a blend cutline from the name OGR support datasource.</dd>
170<dt> <b>-cl</b> <em>layername</em>:</dt><dd>Select the named layer from the
171cutline datasource.</dd>
172<dt> <b>-cwhere</b> <em>expression</em>:</dt><dd>Restrict desired cutline features based on attribute query.</dd>
173<dt> <b>-csql</b> <em>query</em>:</dt><dd>Select cutline features using an SQL query instead of from a layer with -cl.</dd>
174<dt> <b>-cblend</b> <em>distance</em>:</dt><dd>Set a blend distance to use to blend over cutlines (in pixels).</dd>
175<dt> <b>-crop_to_cutline</b>:</dt><dd>(GDAL >= 1.8.0) Crop the extent of the target dataset to the extent of the cutline.</dd>
176<dt> <b>-overwrite</b>:</dt><dd>(GDAL >= 1.8.0) Overwrite the target dataset if it already exists.</dd>
177
178<dt> <em>srcfile</em>:</dt><dd> The source file name(s). </dd>
179<dt> <em>dstfile</em>:</dt><dd> The destination file name. </dd>
180</dl>
181
182Mosaicing into an existing output file is supported if the output file
183already exists. The spatial extent of the existing file will not
184be modified to accomodate new data, so you may have to remove it in that case, or
185use the -overwrite option.
186
187Polygon cutlines may be used as a mask to restrict the area of the destination file
188that may be updated, including blending.  If the OGR layer containing the cutline
189features has no explicit SRS, the cutline features must be in the georeferenced
190units of the destination file. When outputing to a not yet existing target dataset,
191its extent will be the one of the original raster unless -te or -crop_to_cutline are
192specified.
193
194<p>
195\section gdalwarp_example EXAMPLE
196
197For instance, an eight bit spot scene stored in GeoTIFF with
198control points mapping the corners to lat/long could be warped to a UTM
199projection with a command like this:<p>
200
201\verbatim
202gdalwarp -t_srs '+proj=utm +zone=11 +datum=WGS84' raw_spot.tif utm11.tif
203\endverbatim
204
205For instance, the second channel of an ASTER image stored in HDF with
206control points mapping the corners to lat/long could be warped to a UTM
207projection with a command like this:<p>
208
209\verbatim
210gdalwarp HDF4_SDS:ASTER_L1B:"pg-PR1B0000-2002031402_100_001":2 pg-PR1B0000-2002031402_100_001_2.tif
211\endverbatim
212
213\if man
214\section gdalwarp_author AUTHORS
215Frank Warmerdam <warmerdam@pobox.com>, Silke Reimer <silke@intevation.de>
216\endif
217*/
218
219/************************************************************************/
220/*                               GDALExit()                             */
221/*  This function exits and cleans up GDAL and OGR resources            */
222/*  Perhaps it should be added to C api and used in all apps?           */
223/************************************************************************/
224
225static int GDALExit( int nCode )
226{
227  const char  *pszDebug = CPLGetConfigOption("CPL_DEBUG",NULL);
228  if( pszDebug && (EQUAL(pszDebug,"ON") || EQUAL(pszDebug,"") ) )
229  { 
230    GDALDumpOpenDatasets( stderr );
231    CPLDumpSharedList( NULL );
232  }
233
234  GDALDestroyDriverManager();
235
236#ifdef OGR_ENABLED
237  OGRCleanupAll();
238#endif
239
240  return SERVICE_FAILED;
241}
242
243/************************************************************************/
244/*                               Usage()                                */
245/************************************************************************/
246
247static void Usage()
248
249{
250    printf( 
251        "Usage: gdalwarp [--help-general] [--formats]\n"
252        "    [-s_srs srs_def] [-t_srs srs_def] [-to \"NAME=VALUE\"]\n"
253        "    [-order n | -tps | -rpc | -geoloc] [-et err_threshold]\n"
254        "    [-refine_gcps tolerance [minimum_gcps]]\n"
255        "    [-te xmin ymin xmax ymax] [-tr xres yres] [-tap] [-ts width height]\n"
256        "    [-wo \"NAME=VALUE\"] [-ot Byte/Int16/...] [-wt Byte/Int16]\n"
257        "    [-srcnodata \"value [value...]\"] [-dstnodata \"value [value...]\"] -dstalpha\n" 
258        "    [-r resampling_method] [-wm memory_in_mb] [-multi] [-q]\n"
259        "    [-cutline datasource] [-cl layer] [-cwhere expression]\n"
260        "    [-csql statement] [-cblend dist_in_pixels] [-crop_to_cutline]\n"
261        "    [-of format] [-co \"NAME=VALUE\"]* [-overwrite]\n"
262        "    srcfile* dstfile\n"
263        "\n"
264        "Available resampling methods:\n"
265        "    near (default), bilinear, cubic, cubicspline, lanczos.\n" );
266    GDALExit( 1 );
267}
268
269/************************************************************************/
270/*                             SanitizeSRS                              */
271/************************************************************************/
272
273char *SanitizeSRS( const char *pszUserInput )
274
275{
276    OGRSpatialReferenceH hSRS;
277    char *pszResult = NULL;
278
279    CPLErrorReset();
280   
281    hSRS = OSRNewSpatialReference( NULL );
282    if( OSRSetFromUserInput( hSRS, pszUserInput ) == OGRERR_NONE )
283        OSRExportToWkt( hSRS, &pszResult );
284    else
285    {
286        CPLError( CE_Failure, CPLE_AppDefined,
287                  "Translating source or target SRS failed:\n%s",
288                  pszUserInput );
289        GDALExit( 1 );
290    }
291   
292    OSRDestroySpatialReference( hSRS );
293
294    return pszResult;
295}
296
297/************************************************************************/
298/*                                main()                                */
299/************************************************************************/
300
301#ifdef WIN32
302__declspec(dllexport)
303#endif
304int Gdal_Warp( maps*& conf,maps*& inputs,maps*& outputs )
305
306{
307    GDALDatasetH        hDstDS;
308    const char         *pszFormat = "GTiff";
309    int bFormatExplicitelySet = FALSE;
310    char              **papszSrcFiles = NULL;
311    char               *pszDstFilename = NULL;
312    int                 bCreateOutput = FALSE, i;
313    void               *hTransformArg, *hGenImgProjArg=NULL, *hApproxArg=NULL;
314    char               **papszWarpOptions = NULL;
315    double             dfErrorThreshold = 0.125;
316    double             dfWarpMemoryLimit = 0.0;
317    GDALTransformerFunc pfnTransformer = NULL;
318    char                **papszCreateOptions = NULL;
319    GDALDataType        eOutputType = GDT_Unknown, eWorkingType = GDT_Unknown; 
320    GDALResampleAlg     eResampleAlg = GRA_NearestNeighbour;
321    const char          *pszSrcNodata = NULL;
322    const char          *pszDstNodata = NULL;
323    int                 bMulti = FALSE;
324    char                **papszTO = NULL;
325    char                *pszCutlineDSName = NULL;
326    char                *pszCLayer = NULL, *pszCWHERE = NULL, *pszCSQL = NULL;
327    void                *hCutline = NULL;
328    int                  bHasGotErr = FALSE;
329    int                  bCropToCutline = FALSE;
330    int                  bOverwrite = FALSE;
331
332
333    /* Check that we are running against at least GDAL 1.6 */
334    /* Note to developers : if we use newer API, please change the requirement */
335    if (atoi(GDALVersionInfo("VERSION_NUM")) < 1600)
336    {
337        fprintf(stderr, "At least, GDAL >= 1.6.0 is required for this version of warp_service.zo, "
338                "which was compiled against GDAL %s\n", GDAL_RELEASE_NAME);
339        GDALExit(1);
340    }
341
342/* -------------------------------------------------------------------- */
343/*      Register standard GDAL drivers, and process generic GDAL        */
344/*      command options.                                                */
345/* -------------------------------------------------------------------- */
346    GDALAllRegister();
347
348/* -------------------------------------------------------------------- */
349/*      Parse arguments.                                                */
350/* -------------------------------------------------------------------- */
351    char *dataPath;
352    map* tmpMap=getMapFromMaps(conf,"main","dataPath");
353    if(tmpMap!=NULL)
354      dataPath=strdup(tmpMap->value);
355    tmpMap=NULL;
356
357    char *tempPath;
358    tmpMap=getMapFromMaps(conf,"main","tmpPath");
359    if(tmpMap!=NULL){
360      tempPath=strdup(tmpMap->value);
361    }
362
363    tmpMap=NULL;
364    tmpMap=getMapFromMaps(inputs,"InputDSN","value");
365    if(tmpMap!=NULL){
366      map* length=getMapFromMaps(inputs,"InputDSN","length");
367      if(length!=NULL){
368        int i=0;
369        int len=atoi(length->value);
370        maps* currentMaps=getMaps(inputs,"InputDSN");
371        for(i=0;i<len;i++){
372          tmpMap=getMapArray(currentMaps->content,"value",i);
373          papszSrcFiles[i]=(char*)CPLMalloc(sizeof(char)*(strlen(dataPath)+strlen(tmpMap->value)+4));
374          sprintf((char*)papszSrcFiles[i],"%s/%s",dataPath,tmpMap->value);
375        }
376      }else{
377        char *tmp=(char*)CPLMalloc(sizeof(char)*(strlen(dataPath)+strlen(tmpMap->value)+4));
378        sprintf(tmp,"%s/%s",dataPath,tmpMap->value);
379        papszSrcFiles=CSLAddString( papszSrcFiles, tmp );
380        free(tmp);
381      }
382    }
383
384    tmpMap=NULL;
385    tmpMap=getMapFromMaps(inputs,"OutputDSN","value");
386    if(tmpMap!=NULL){
387      pszDstFilename=(char*)CPLMalloc(sizeof(char)*(strlen(tempPath)+strlen(tmpMap->value)+4));
388      char *ext=new char[4];
389      ext="tif";
390      if(strncasecmp(pszFormat,"AAIGRID",7)==0)
391        ext="csv";
392      else 
393        if(strncasecmp(pszFormat,"PNG",3)==0)
394          ext="png";
395        else
396          if(strncasecmp(pszFormat,"GIF",3)==0)
397            ext="gif";
398          else
399            if(strncasecmp(pszFormat,"JPEG",4)==0)
400              ext="jpg";
401      sprintf(pszDstFilename,"%s/%s.%s",tempPath,tmpMap->value,ext);
402    }
403
404
405    tmpMap=NULL;
406    tmpMap=getMapFromMaps(inputs,"r","value");
407    if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
408      if ( EQUAL(tmpMap->value, "near") )
409        eResampleAlg = GRA_NearestNeighbour;
410      else if ( EQUAL(tmpMap->value, "bilinear") )
411        eResampleAlg = GRA_Bilinear;
412      else if ( EQUAL(tmpMap->value, "cubic") )
413        eResampleAlg = GRA_Cubic;
414      else if ( EQUAL(tmpMap->value, "cubicspline") )
415        eResampleAlg = GRA_CubicSpline;
416      else if ( EQUAL(tmpMap->value, "lanczos") )
417        eResampleAlg = GRA_Lanczos;
418      else
419        {
420          char tmp[1024];
421          sprintf(tmp,"Unknown resampling method: \"%s\".", tmpMap->value );
422          setMapInMaps(conf,"lenv","message",tmp);
423          return SERVICE_FAILED;
424        }
425    }
426
427    tmpMap=NULL;
428    tmpMap=getMapFromMaps(inputs,"t_srs","value");
429    if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
430      char *pszSRS = SanitizeSRS(tmpMap->value);
431      papszTO = CSLSetNameValue( papszTO, "DST_SRS", pszSRS );
432      CPLFree( pszSRS );
433    }
434
435    tmpMap=getMapFromMaps(inputs,"s_srs","value");
436    if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
437      char *pszSRS = SanitizeSRS(tmpMap->value);
438      papszTO = CSLSetNameValue( papszTO, "SRC_SRS", pszSRS );
439      CPLFree( pszSRS );
440    }
441
442    tmpMap=getMapFromMaps(inputs,"order","value");
443    if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
444      const char* pszMethod = CSLFetchNameValue(papszTO, "METHOD");
445      if (pszMethod)
446        fprintf(stderr, "Warning: only one METHOD can be used. Method %s is already defined\n",
447                pszMethod);
448      papszTO = CSLSetNameValue( papszTO, "MAX_GCP_ORDER", tmpMap->value );
449    }
450   
451    tmpMap=getMapFromMaps(inputs,"co","value");
452    if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
453      papszCreateOptions = CSLAddString( papszCreateOptions, tmpMap->value );
454      bCreateOutput = TRUE;
455    }   
456
457    tmpMap=getMapFromMaps(inputs,"wo","value");
458    if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
459      papszWarpOptions = CSLAddString( papszWarpOptions, tmpMap->value );
460    }   
461
462    tmpMap=getMapFromMaps(inputs,"multi","value");
463    if(tmpMap!=NULL && strncasecmp(tmpMap->value,"TRUE",4)==0 ){
464      bMulti = TRUE;
465    }   
466
467    tmpMap=getMapFromMaps(inputs,"dstalpha","value");
468    if(tmpMap!=NULL && strncasecmp(tmpMap->value,"TRUE",4)==0 ){
469      bEnableDstAlpha = TRUE;
470    }
471
472    tmpMap=getMapFromMaps(inputs,"srcalpha","value");
473    if(tmpMap!=NULL && strncasecmp(tmpMap->value,"TRUE",4)==0 ){
474      bEnableSrcAlpha = TRUE;
475    }
476
477    tmpMap=getMapFromMaps(inputs,"tap","value");
478    if(tmpMap!=NULL && strncasecmp(tmpMap->value,"TRUE",4)==0 ){
479      bTargetAlignedPixels = TRUE;
480    }
481
482
483    tmpMap=getMapFromMaps(inputs,"crop_to_cutline","value");
484    if(tmpMap!=NULL && strncasecmp(tmpMap->value,"TRUE",4)==0 ){
485      bCropToCutline = TRUE;
486      bCreateOutput = TRUE;
487    }
488
489
490    tmpMap=getMapFromMaps(inputs,"tps","value");
491    if(tmpMap!=NULL && strncasecmp(tmpMap->value,"TRUE",4)==0 ){
492      papszTO = CSLSetNameValue( papszTO, "METHOD", "GCP_TPS" );
493    }
494
495    tmpMap=getMapFromMaps(inputs,"rpc","value");
496    if(tmpMap!=NULL && strncasecmp(tmpMap->value,"TRUE",4)==0 ){
497      papszTO = CSLSetNameValue( papszTO, "METHOD", "RPC" );
498    }
499
500    tmpMap=getMapFromMaps(inputs,"geoloc","value");
501    if(tmpMap!=NULL && strncasecmp(tmpMap->value,"TRUE",4)==0 ){
502      papszTO = CSLSetNameValue( papszTO, "METHOD", "GEOLOC_ARRAY" );
503    }
504
505    tmpMap=getMapFromMaps(inputs,"wt","value");
506    if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
507      int       iType;
508     
509      for( iType = 1; iType < GDT_TypeCount; iType++ )
510        {
511          if( GDALGetDataTypeName((GDALDataType)iType) != NULL
512              && EQUAL(GDALGetDataTypeName((GDALDataType)iType),
513                       tmpMap->value) )
514            {
515              eWorkingType = (GDALDataType) iType;
516            }
517        }
518     
519      if( eWorkingType == GDT_Unknown )
520        {
521          char tmp[1024];
522          sprintf(tmp,"Unknown output pixel type: %s\n", tmpMap->value );
523          setMapInMaps(conf,"lenv","message",tmp);
524          return SERVICE_FAILED;
525        }
526    }
527
528    tmpMap=getMapFromMaps(inputs,"srcnodata","value");
529    if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
530      pszSrcNodata = strdup(tmpMap->value);
531    }
532
533    tmpMap=getMapFromMaps(inputs,"dstnodata","value");
534    if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
535      pszDstNodata = strdup(tmpMap->value);
536    }
537
538
539    tmpMap=getMapFromMaps(inputs,"tr","value");
540    if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
541      char *tmpV=strdup(tmpMap->value);
542      char *res=strtok(tmpV,",");
543      int i=0;
544      while(res!=NULL){
545        if(i==0)
546          dfXRes = CPLAtofM(res);
547        else
548          dfYRes = fabs(CPLAtofM(res));
549        res=strtok(NULL,",");
550        i++;
551      }
552      if( dfXRes == 0 || dfYRes == 0 )
553        {
554          setMapInMaps(conf,"lenv","message","Wrong value for TR parameters");
555        }
556      bCreateOutput = TRUE;
557    }
558
559    tmpMap=getMapFromMaps(inputs,"cblend","value");
560    if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
561      papszWarpOptions = 
562        CSLSetNameValue( papszWarpOptions, 
563                         "CUTLINE_BLEND_DIST", tmpMap->value );
564    }
565   
566    tmpMap=getMapFromMaps(inputs,"ot","value");
567    if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
568      int       iType;
569     
570      for( iType = 1; iType < GDT_TypeCount; iType++ )
571        {
572          if( GDALGetDataTypeName((GDALDataType)iType) != NULL
573              && EQUAL(GDALGetDataTypeName((GDALDataType)iType),
574                       tmpMap->value) )
575            {
576              eOutputType = (GDALDataType) iType;
577            }
578        }
579
580      if( eOutputType == GDT_Unknown )
581        {
582          char tmp[128];
583          sprintf( tmp,"Unknown output pixel type: %s\n", tmpMap->value );
584          setMapInMaps(conf,"lenv","message",tmp);
585          return SERVICE_FAILED;
586        }
587      bCreateOutput = TRUE;     
588    }
589
590
591    tmpMap=getMapFromMaps(inputs,"ts","value");
592    if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
593      char *tmpV=strdup(tmpMap->value);
594      char *res=strtok(tmpV,",");
595      int i=0;
596      while(res!=NULL){
597        if(i==0)
598          nForcePixels = atoi(res);       
599        else
600          nForceLines = atoi(res);       
601        res=strtok(NULL,",");
602        i++;
603      }
604      free(tmpV);
605      bCreateOutput = TRUE;     
606    }
607
608    tmpMap=getMapFromMaps(inputs,"te","value");
609    if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
610      char *tmpV=strdup(tmpMap->value);
611      char *res=strtok(tmpV,",");
612      int i=0;
613      while(res!=NULL){
614        switch(i){
615        case 0:
616          dfMinX = CPLAtofM(res);
617          break;
618        case 1:
619          dfMinY = CPLAtofM(res);
620          break;
621        case 2:
622          dfMaxX = CPLAtofM(res);
623          break;
624        case 3:
625          dfMaxY = CPLAtofM(res);
626          break;
627        }
628        if(i==0)
629          nForcePixels = atoi(res);       
630        else
631          nForceLines = atoi(res);       
632        res=strtok(NULL,",");
633        i++;
634      }
635      free(tmpV);
636      bCreateOutput = TRUE;
637    }
638
639    tmpMap=getMapFromMaps(inputs,"te","value");
640    if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
641      pszCutlineDSName = strdup(tmpMap->value);
642    }
643
644    tmpMap=getMapFromMaps(inputs,"cwhere","value");
645    if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
646      pszCWHERE = strdup(tmpMap->value);
647    }
648
649    tmpMap=getMapFromMaps(inputs,"cl","value");
650    if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
651      pszCLayer = strdup(tmpMap->value);
652    }
653
654    tmpMap=getMapFromMaps(inputs,"csql","value");
655    if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
656      pszCSQL = strdup(tmpMap->value);
657    }
658   
659    tmpMap=getMapFromMaps(inputs,"of","value");
660    if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
661      pszFormat = strdup(tmpMap->value);
662      bFormatExplicitelySet = TRUE;
663      bCreateOutput = TRUE;
664      if( EQUAL(pszFormat,"VRT") )
665        bVRT = TRUE;
666    }
667
668    tmpMap=getMapFromMaps(inputs,"to","value");
669    if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
670      papszTO = CSLAddString( papszTO, tmpMap->value );
671    }
672
673    tmpMap=getMapFromMaps(inputs,"et","value");
674    if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
675      dfErrorThreshold = CPLAtofM(tmpMap->value);
676    }
677   
678
679    tmpMap=getMapFromMaps(inputs,"refine_gpcs","value");
680    if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
681      char *tmpV=strdup(tmpMap->value);
682      char *res=strtok(tmpV,",");
683      int i=0;
684      while(res!=NULL){
685        if(i==0){
686          papszTO = CSLSetNameValue( papszTO, "REFINE_TOLERANCE", res );
687          if(atof(res) < 0){
688            setMapInMaps(conf,"lenv","message","The tolerance for -refine_gcps may not be negative\n");
689            return SERVICE_FAILED;
690          }
691        }else if(i==1){
692          if(atoi(res) >= 0 && isdigit(res[0]))
693            papszTO = CSLSetNameValue( papszTO, "REFINE_MINIMUM_GCPS", res );
694          else
695            papszTO = CSLSetNameValue( papszTO, "REFINE_MINIMUM_GCPS", "-1" );
696        }
697        res=strtok(NULL,",");
698        i++;
699      }
700    }
701
702    bQuiet = TRUE;
703   
704/* -------------------------------------------------------------------- */
705/*      Check that incompatible options are not used                    */
706/* -------------------------------------------------------------------- */
707
708    if ((nForcePixels != 0 || nForceLines != 0) && 
709        (dfXRes != 0 && dfYRes != 0))
710    {
711      setMapInMaps(conf,"lenv","message","tr and ts options cannot be used at the same time\n");
712      return SERVICE_FAILED;
713    }
714   
715    if (bTargetAlignedPixels && dfXRes == 0 && dfYRes == 0)
716    {
717      setMapInMaps(conf,"lenv","message","tap option cannot be used without using tr");
718      return SERVICE_FAILED;
719    }
720
721
722/* -------------------------------------------------------------------- */
723/*      Does the output dataset already exist?                          */
724/* -------------------------------------------------------------------- */
725
726    /* FIXME ? source filename=target filename and -overwrite is definitely */
727    /* an error. But I can't imagine of a valid case (without -overwrite), */
728    /* where it would make sense. In doubt, let's keep that dubious possibility... */
729    if ( CSLCount(papszSrcFiles) == 1 &&
730         strcmp(papszSrcFiles[0], pszDstFilename) == 0 && bOverwrite)
731    {
732      setMapInMaps(conf,"lenv","message","Source and destination datasets must be different.");
733      return GDALExit(1);
734    }
735
736    CPLPushErrorHandler( CPLQuietErrorHandler );
737    hDstDS = GDALOpen( pszDstFilename, GA_Update );
738    CPLPopErrorHandler();
739
740    if( hDstDS != NULL && bOverwrite )
741    {
742        GDALClose(hDstDS);
743        hDstDS = NULL;
744    }
745
746    if( hDstDS != NULL && bCreateOutput )
747    {
748      char tmp[024];
749      sprintf( tmp, 
750               "Output dataset %s exists,\n"
751               "but some commandline options were provided indicating a new dataset\n"
752               "should be created.  Please delete existing dataset and run again.\n",
753               pszDstFilename );
754      setMapInMaps(conf,"lenv","message",tmp);
755      return GDALExit(1);
756    }
757
758    /* Avoid overwriting an existing destination file that cannot be opened in */
759    /* update mode with a new GTiff file */
760    if ( hDstDS == NULL && !bOverwrite )
761    {
762        CPLPushErrorHandler( CPLQuietErrorHandler );
763        hDstDS = GDALOpen( pszDstFilename, GA_ReadOnly );
764        CPLPopErrorHandler();
765       
766        if (hDstDS)
767        {
768          char tmp[1024];
769          sprintf( tmp, 
770                   "Output dataset %s exists, but cannot be opened in update mode\n",
771                   pszDstFilename );
772          GDALClose(hDstDS);
773          setMapInMaps(conf,"lenv","message",tmp);
774          return GDALExit( 1 );
775        }
776    }
777
778/* -------------------------------------------------------------------- */
779/*      If we have a cutline datasource read it and attach it in the    */
780/*      warp options.                                                   */
781/* -------------------------------------------------------------------- */
782    if( pszCutlineDSName != NULL )
783    {
784        LoadCutline( pszCutlineDSName, pszCLayer, pszCWHERE, pszCSQL,
785                     &hCutline );
786    }
787
788#ifdef OGR_ENABLED
789    if ( bCropToCutline && hCutline != NULL )
790    {
791     
792     
793        OGRGeometryH hCutlineGeom = OGR_G_Clone( (OGRGeometryH) hCutline );
794        OGRSpatialReferenceH hCutlineSRS = OGR_G_GetSpatialReference( hCutlineGeom );
795        const char *pszThisTargetSRS = CSLFetchNameValue( papszTO, "DST_SRS" );
796        OGRCoordinateTransformationH hCT = NULL;
797        if (hCutlineSRS == NULL)
798        {
799            /* We suppose it is in target coordinates */
800        }
801        else if (pszThisTargetSRS != NULL)
802        {
803            OGRSpatialReferenceH hTargetSRS = OSRNewSpatialReference(NULL);
804            if( OSRImportFromWkt( hTargetSRS, (char **)&pszThisTargetSRS ) != CE_None )
805            {
806              setMapInMaps(conf,"lenv","message","Cannot compute bounding box of cutline.\n");
807                return GDALExit(1);
808            }
809
810            hCT = OCTNewCoordinateTransformation(hCutlineSRS, hTargetSRS);
811
812            OSRDestroySpatialReference(hTargetSRS);
813        }
814        else if (pszThisTargetSRS == NULL)
815        {
816            if (papszSrcFiles[0] != NULL)
817            {
818                GDALDatasetH hSrcDS = GDALOpen(papszSrcFiles[0], GA_ReadOnly);
819                if (hSrcDS == NULL)
820                {
821                  setMapInMaps(conf,"lenv","message","Cannot compute bounding box of cutline.");
822                  return GDALExit(1);
823                }
824
825                OGRSpatialReferenceH  hRasterSRS = NULL;
826                const char *pszProjection = NULL;
827
828                if( GDALGetProjectionRef( hSrcDS ) != NULL
829                    && strlen(GDALGetProjectionRef( hSrcDS )) > 0 )
830                    pszProjection = GDALGetProjectionRef( hSrcDS );
831                else if( GDALGetGCPProjection( hSrcDS ) != NULL )
832                    pszProjection = GDALGetGCPProjection( hSrcDS );
833
834                if( pszProjection == NULL )
835                {
836                  setMapInMaps(conf,"lenv","message","Cannot compute bounding box of cutline.");
837                  return GDALExit(1);
838                }
839
840                hRasterSRS = OSRNewSpatialReference(NULL);
841                if( OSRImportFromWkt( hRasterSRS, (char **)&pszProjection ) != CE_None )
842                {
843                  setMapInMaps(conf,"lenv","message","Cannot compute bounding box of cutline.");
844                  return GDALExit(1);
845                }
846
847                hCT = OCTNewCoordinateTransformation(hCutlineSRS, hRasterSRS);
848
849                OSRDestroySpatialReference(hRasterSRS);
850
851                GDALClose(hSrcDS);
852            }
853            else
854            {
855              setMapInMaps(conf,"levn","message", "Cannot compute bounding box of cutline.\n");
856              return GDALExit(1);
857            }
858        }
859
860        if (hCT)
861        {
862            OGR_G_Transform( hCutlineGeom, hCT );
863
864            OCTDestroyCoordinateTransformation(hCT);
865        }
866
867        OGREnvelope sEnvelope;
868        OGR_G_GetEnvelope(hCutlineGeom, &sEnvelope);
869
870        dfMinX = sEnvelope.MinX;
871        dfMinY = sEnvelope.MinY;
872        dfMaxX = sEnvelope.MaxX;
873        dfMaxY = sEnvelope.MaxY;
874       
875        OGR_G_DestroyGeometry(hCutlineGeom);
876    }
877#endif
878   
879/* -------------------------------------------------------------------- */
880/*      If not, we need to create it.                                   */
881/* -------------------------------------------------------------------- */
882    int   bInitDestSetForFirst = FALSE;
883
884    void* hUniqueTransformArg = NULL;
885    GDALDatasetH hUniqueSrcDS = NULL;
886
887    if( hDstDS == NULL )
888    {
889
890      hDstDS = GDALWarpCreateOutput( conf,papszSrcFiles, pszDstFilename,pszFormat,
891                                       papszTO, &papszCreateOptions, 
892                                       eOutputType, &hUniqueTransformArg,
893                                       &hUniqueSrcDS);
894        bCreateOutput = TRUE;
895
896
897        if( CSLFetchNameValue( papszWarpOptions, "INIT_DEST" ) == NULL 
898            && pszDstNodata == NULL )
899        {
900            papszWarpOptions = CSLSetNameValue(papszWarpOptions,
901                                               "INIT_DEST", "0");
902            bInitDestSetForFirst = TRUE;
903        }
904        else if( CSLFetchNameValue( papszWarpOptions, "INIT_DEST" ) == NULL )
905        {
906            papszWarpOptions = CSLSetNameValue(papszWarpOptions,
907                                               "INIT_DEST", "NO_DATA" );
908            bInitDestSetForFirst = TRUE;
909        }
910
911        CSLDestroy( papszCreateOptions );
912        papszCreateOptions = NULL;
913    }
914
915 
916    if( hDstDS == NULL )
917        GDALExit( 1 );
918
919/* -------------------------------------------------------------------- */
920/*      Loop over all source files, processing each in turn.            */
921/* -------------------------------------------------------------------- */
922    int iSrc;
923
924    for( iSrc = 0; papszSrcFiles[iSrc] != NULL; iSrc++ )
925    {
926        GDALDatasetH hSrcDS;
927       
928/* -------------------------------------------------------------------- */
929/*      Open this file.                                                 */
930/* -------------------------------------------------------------------- */
931        if (hUniqueSrcDS)
932            hSrcDS = hUniqueSrcDS;
933        else
934            hSrcDS = GDALOpen( papszSrcFiles[iSrc], GA_ReadOnly );
935   
936        if( hSrcDS == NULL ){
937          setMapInMaps(conf,"lenv","message",CPLGetLastErrorMsg());
938          return GDALExit( 2 );
939        }
940
941/* -------------------------------------------------------------------- */
942/*      Check that there's at least one raster band                     */
943/* -------------------------------------------------------------------- */
944        if ( GDALGetRasterCount(hSrcDS) == 0 )
945        {   
946          char tmp[1024];
947          sprintf(tmp, "Input file %s has no raster bands.\n", papszSrcFiles[iSrc] );
948          setMapInMaps(conf,"lenv","message",tmp);
949          return GDALExit( 1 );
950        }
951       
952        if( !bQuiet )
953            printf( "Processing input file %s.\n", papszSrcFiles[iSrc] );
954
955/* -------------------------------------------------------------------- */
956/*      Warns if the file has a color table and something more          */
957/*      complicated than nearest neighbour resampling is asked          */
958/* -------------------------------------------------------------------- */
959 
960        if ( eResampleAlg != GRA_NearestNeighbour &&
961             GDALGetRasterColorTable(GDALGetRasterBand(hSrcDS, 1)) != NULL)
962        {
963                fprintf( stderr, "Warning: Input file %s has a color table, which will likely lead to "
964                        "bad results when using a resampling method other than "
965                        "nearest neighbour. Converting the dataset prior to 24/32 bit "
966                        "is advised.\n", papszSrcFiles[iSrc] );
967        }
968
969/* -------------------------------------------------------------------- */
970/*      Do we have a source alpha band?                                 */
971/* -------------------------------------------------------------------- */
972        if( GDALGetRasterColorInterpretation( 
973                GDALGetRasterBand(hSrcDS,GDALGetRasterCount(hSrcDS)) ) 
974            == GCI_AlphaBand
975            && !bEnableSrcAlpha )
976        {
977            bEnableSrcAlpha = TRUE;
978            if( !bQuiet )
979                printf( "Using band %d of source image as alpha.\n", 
980                        GDALGetRasterCount(hSrcDS) );
981        }
982
983/* -------------------------------------------------------------------- */
984/*      Create a transformation object from the source to               */
985/*      destination coordinate system.                                  */
986/* -------------------------------------------------------------------- */
987        if (hUniqueTransformArg)
988            hTransformArg = hGenImgProjArg = hUniqueTransformArg;
989        else
990            hTransformArg = hGenImgProjArg =
991                GDALCreateGenImgProjTransformer2( hSrcDS, hDstDS, papszTO );
992
993        if( hTransformArg == NULL ){
994          setMapInMaps(conf,"lenv","message",CPLGetLastErrorMsg());
995          return GDALExit( 1 );
996        }
997       
998        pfnTransformer = GDALGenImgProjTransform;
999
1000/* -------------------------------------------------------------------- */
1001/*      Warp the transformer with a linear approximator unless the      */
1002/*      acceptable error is zero.                                       */
1003/* -------------------------------------------------------------------- */
1004        if( dfErrorThreshold != 0.0 )
1005        {
1006            hTransformArg = hApproxArg = 
1007                GDALCreateApproxTransformer( GDALGenImgProjTransform, 
1008                                             hGenImgProjArg, dfErrorThreshold);
1009            pfnTransformer = GDALApproxTransform;
1010        }
1011
1012/* -------------------------------------------------------------------- */
1013/*      Clear temporary INIT_DEST settings after the first image.       */
1014/* -------------------------------------------------------------------- */
1015        if( bInitDestSetForFirst && iSrc == 1 )
1016            papszWarpOptions = CSLSetNameValue( papszWarpOptions, 
1017                                                "INIT_DEST", NULL );
1018
1019/* -------------------------------------------------------------------- */
1020/*      Setup warp options.                                             */
1021/* -------------------------------------------------------------------- */
1022        GDALWarpOptions *psWO = GDALCreateWarpOptions();
1023
1024        psWO->papszWarpOptions = CSLDuplicate(papszWarpOptions);
1025        psWO->eWorkingDataType = eWorkingType;
1026        psWO->eResampleAlg = eResampleAlg;
1027
1028        psWO->hSrcDS = hSrcDS;
1029        psWO->hDstDS = hDstDS;
1030
1031        psWO->pfnTransformer = pfnTransformer;
1032        psWO->pTransformerArg = hTransformArg;
1033
1034        if( !bQuiet )
1035            psWO->pfnProgress = GDALTermProgress;
1036
1037        if( dfWarpMemoryLimit != 0.0 )
1038            psWO->dfWarpMemoryLimit = dfWarpMemoryLimit;
1039
1040/* -------------------------------------------------------------------- */
1041/*      Setup band mapping.                                             */
1042/* -------------------------------------------------------------------- */
1043        if( bEnableSrcAlpha )
1044            psWO->nBandCount = GDALGetRasterCount(hSrcDS) - 1;
1045        else
1046            psWO->nBandCount = GDALGetRasterCount(hSrcDS);
1047
1048        psWO->panSrcBands = (int *) CPLMalloc(psWO->nBandCount*sizeof(int));
1049        psWO->panDstBands = (int *) CPLMalloc(psWO->nBandCount*sizeof(int));
1050
1051        for( i = 0; i < psWO->nBandCount; i++ )
1052        {
1053            psWO->panSrcBands[i] = i+1;
1054            psWO->panDstBands[i] = i+1;
1055        }
1056
1057/* -------------------------------------------------------------------- */
1058/*      Setup alpha bands used if any.                                  */
1059/* -------------------------------------------------------------------- */
1060        if( bEnableSrcAlpha )
1061            psWO->nSrcAlphaBand = GDALGetRasterCount(hSrcDS);
1062
1063        if( !bEnableDstAlpha
1064            && GDALGetRasterCount(hDstDS) == psWO->nBandCount+1 
1065            && GDALGetRasterColorInterpretation( 
1066                GDALGetRasterBand(hDstDS,GDALGetRasterCount(hDstDS))) 
1067            == GCI_AlphaBand )
1068        {
1069            if( !bQuiet )
1070                printf( "Using band %d of destination image as alpha.\n", 
1071                        GDALGetRasterCount(hDstDS) );
1072               
1073            bEnableDstAlpha = TRUE;
1074        }
1075
1076        if( bEnableDstAlpha )
1077            psWO->nDstAlphaBand = GDALGetRasterCount(hDstDS);
1078
1079/* -------------------------------------------------------------------- */
1080/*      Setup NODATA options.                                           */
1081/* -------------------------------------------------------------------- */
1082        if( pszSrcNodata != NULL && !EQUALN(pszSrcNodata,"n",1) )
1083        {
1084            char **papszTokens = CSLTokenizeString( pszSrcNodata );
1085            int  nTokenCount = CSLCount(papszTokens);
1086
1087            psWO->padfSrcNoDataReal = (double *) 
1088                CPLMalloc(psWO->nBandCount*sizeof(double));
1089            psWO->padfSrcNoDataImag = (double *) 
1090                CPLMalloc(psWO->nBandCount*sizeof(double));
1091
1092            for( i = 0; i < psWO->nBandCount; i++ )
1093            {
1094                if( i < nTokenCount )
1095                {
1096                    CPLStringToComplex( papszTokens[i], 
1097                                        psWO->padfSrcNoDataReal + i,
1098                                        psWO->padfSrcNoDataImag + i );
1099                }
1100                else
1101                {
1102                    psWO->padfSrcNoDataReal[i] = psWO->padfSrcNoDataReal[i-1];
1103                    psWO->padfSrcNoDataImag[i] = psWO->padfSrcNoDataImag[i-1];
1104                }
1105            }
1106
1107            CSLDestroy( papszTokens );
1108
1109            psWO->papszWarpOptions = CSLSetNameValue(psWO->papszWarpOptions,
1110                                               "UNIFIED_SRC_NODATA", "YES" );
1111        }
1112
1113/* -------------------------------------------------------------------- */
1114/*      If -srcnodata was not specified, but the data has nodata        */
1115/*      values, use them.                                               */
1116/* -------------------------------------------------------------------- */
1117        if( pszSrcNodata == NULL )
1118        {
1119            int bHaveNodata = FALSE;
1120            double dfReal = 0.0;
1121
1122            for( i = 0; !bHaveNodata && i < psWO->nBandCount; i++ )
1123            {
1124                GDALRasterBandH hBand = GDALGetRasterBand( hSrcDS, i+1 );
1125                dfReal = GDALGetRasterNoDataValue( hBand, &bHaveNodata );
1126            }
1127
1128            if( bHaveNodata )
1129            {
1130                if( !bQuiet )
1131                {
1132                    if (CPLIsNan(dfReal))
1133                        printf( "Using internal nodata values (eg. nan) for image %s.\n",
1134                                papszSrcFiles[iSrc] );
1135                    else
1136                        printf( "Using internal nodata values (eg. %g) for image %s.\n",
1137                                dfReal, papszSrcFiles[iSrc] );
1138                }
1139                psWO->padfSrcNoDataReal = (double *) 
1140                    CPLMalloc(psWO->nBandCount*sizeof(double));
1141                psWO->padfSrcNoDataImag = (double *) 
1142                    CPLMalloc(psWO->nBandCount*sizeof(double));
1143               
1144                for( i = 0; i < psWO->nBandCount; i++ )
1145                {
1146                    GDALRasterBandH hBand = GDALGetRasterBand( hSrcDS, i+1 );
1147
1148                    dfReal = GDALGetRasterNoDataValue( hBand, &bHaveNodata );
1149
1150                    if( bHaveNodata )
1151                    {
1152                        psWO->padfSrcNoDataReal[i] = dfReal;
1153                        psWO->padfSrcNoDataImag[i] = 0.0;
1154                    }
1155                    else
1156                    {
1157                        psWO->padfSrcNoDataReal[i] = -123456.789;
1158                        psWO->padfSrcNoDataImag[i] = 0.0;
1159                    }
1160                }
1161            }
1162        }
1163
1164/* -------------------------------------------------------------------- */
1165/*      If the output dataset was created, and we have a destination    */
1166/*      nodata value, go through marking the bands with the information.*/
1167/* -------------------------------------------------------------------- */
1168        if( pszDstNodata != NULL )
1169        {
1170            char **papszTokens = CSLTokenizeString( pszDstNodata );
1171            int  nTokenCount = CSLCount(papszTokens);
1172
1173            psWO->padfDstNoDataReal = (double *) 
1174                CPLMalloc(psWO->nBandCount*sizeof(double));
1175            psWO->padfDstNoDataImag = (double *) 
1176                CPLMalloc(psWO->nBandCount*sizeof(double));
1177
1178            for( i = 0; i < psWO->nBandCount; i++ )
1179            {
1180                if( i < nTokenCount )
1181                {
1182                    CPLStringToComplex( papszTokens[i], 
1183                                        psWO->padfDstNoDataReal + i,
1184                                        psWO->padfDstNoDataImag + i );
1185                }
1186                else
1187                {
1188                    psWO->padfDstNoDataReal[i] = psWO->padfDstNoDataReal[i-1];
1189                    psWO->padfDstNoDataImag[i] = psWO->padfDstNoDataImag[i-1];
1190                }
1191               
1192                GDALRasterBandH hBand = GDALGetRasterBand( hDstDS, i+1 );
1193                int bClamped = FALSE, bRounded = FALSE;
1194
1195#define CLAMP(val,type,minval,maxval) \
1196    do { if (val < minval) { bClamped = TRUE; val = minval; } \
1197    else if (val > maxval) { bClamped = TRUE; val = maxval; } \
1198    else if (val != (type)val) { bRounded = TRUE; val = (type)(val + 0.5); } } \
1199    while(0)
1200
1201                switch(GDALGetRasterDataType(hBand))
1202                {
1203                    case GDT_Byte:
1204                        CLAMP(psWO->padfDstNoDataReal[i], GByte,
1205                              0.0, 255.0);
1206                        break;
1207                    case GDT_Int16:
1208                        CLAMP(psWO->padfDstNoDataReal[i], GInt16,
1209                              -32768.0, 32767.0);
1210                        break;
1211                    case GDT_UInt16:
1212                        CLAMP(psWO->padfDstNoDataReal[i], GUInt16,
1213                              0.0, 65535.0);
1214                        break;
1215                    case GDT_Int32:
1216                        CLAMP(psWO->padfDstNoDataReal[i], GInt32,
1217                              -2147483648.0, 2147483647.0);
1218                        break;
1219                    case GDT_UInt32:
1220                        CLAMP(psWO->padfDstNoDataReal[i], GUInt32,
1221                              0.0, 4294967295.0);
1222                        break;
1223                    default:
1224                        break;
1225                }
1226                   
1227                if (bClamped)
1228                {
1229                    printf( "for band %d, destination nodata value has been clamped "
1230                           "to %.0f, the original value being out of range.\n",
1231                           i + 1, psWO->padfDstNoDataReal[i]);
1232                }
1233                else if(bRounded)
1234                {
1235                    printf("for band %d, destination nodata value has been rounded "
1236                           "to %.0f, %s being an integer datatype.\n",
1237                           i + 1, psWO->padfDstNoDataReal[i],
1238                           GDALGetDataTypeName(GDALGetRasterDataType(hBand)));
1239                }
1240
1241                if( bCreateOutput )
1242                {
1243                    GDALSetRasterNoDataValue( 
1244                        GDALGetRasterBand( hDstDS, psWO->panDstBands[i] ), 
1245                        psWO->padfDstNoDataReal[i] );
1246                }
1247            }
1248
1249            CSLDestroy( papszTokens );
1250        }
1251
1252/* -------------------------------------------------------------------- */
1253/*      If we have a cutline, transform it into the source              */
1254/*      pixel/line coordinate system and insert into warp options.      */
1255/* -------------------------------------------------------------------- */
1256        if( hCutline != NULL )
1257        {
1258            TransformCutlineToSource( hSrcDS, hCutline, 
1259                                      &(psWO->papszWarpOptions), 
1260                                      papszTO );
1261        }
1262
1263/* -------------------------------------------------------------------- */
1264/*      If we are producing VRT output, then just initialize it with    */
1265/*      the warp options and write out now rather than proceeding       */
1266/*      with the operations.                                            */
1267/* -------------------------------------------------------------------- */
1268        if( bVRT )
1269        {
1270            if( GDALInitializeWarpedVRT( hDstDS, psWO ) != CE_None )
1271                GDALExit( 1 );
1272
1273            GDALClose( hDstDS );
1274            GDALClose( hSrcDS );
1275
1276            /* The warped VRT will clean itself the transformer used */
1277            /* So we have only to destroy the hGenImgProjArg if we */
1278            /* have wrapped it inside the hApproxArg */
1279            if (pfnTransformer == GDALApproxTransform)
1280            {
1281                if( hGenImgProjArg != NULL )
1282                    GDALDestroyGenImgProjTransformer( hGenImgProjArg );
1283            }
1284
1285            GDALDestroyWarpOptions( psWO );
1286
1287            setMapInMaps(outputs,"Result","value",pszDstFilename);
1288
1289            CPLFree( pszDstFilename );
1290            CSLDestroy( papszSrcFiles );
1291            CSLDestroy( papszWarpOptions );
1292            CSLDestroy( papszTO );
1293   
1294            GDALDumpOpenDatasets( stderr );
1295       
1296            GDALDestroyDriverManager();
1297       
1298            return SERVICE_SUCCEEDED;
1299        }
1300
1301/* -------------------------------------------------------------------- */
1302/*      Initialize and execute the warp.                                */
1303/* -------------------------------------------------------------------- */
1304        GDALWarpOperation oWO;
1305
1306        if( oWO.Initialize( psWO ) == CE_None )
1307        {
1308            CPLErr eErr;
1309            if( bMulti )
1310                eErr = oWO.ChunkAndWarpMulti( 0, 0, 
1311                                       GDALGetRasterXSize( hDstDS ),
1312                                       GDALGetRasterYSize( hDstDS ) );
1313            else
1314                eErr = oWO.ChunkAndWarpImage( 0, 0, 
1315                                       GDALGetRasterXSize( hDstDS ),
1316                                       GDALGetRasterYSize( hDstDS ) );
1317            if (eErr != CE_None){
1318                bHasGotErr = TRUE;
1319            }
1320        }
1321
1322/* -------------------------------------------------------------------- */
1323/*      Cleanup                                                         */
1324/* -------------------------------------------------------------------- */
1325        if( hApproxArg != NULL )
1326            GDALDestroyApproxTransformer( hApproxArg );
1327       
1328        if( hGenImgProjArg != NULL )
1329            GDALDestroyGenImgProjTransformer( hGenImgProjArg );
1330       
1331        GDALDestroyWarpOptions( psWO );
1332
1333        GDALClose( hSrcDS );
1334    }
1335
1336/* -------------------------------------------------------------------- */
1337/*      Final Cleanup.                                                  */
1338/* -------------------------------------------------------------------- */
1339    CPLErrorReset();
1340    GDALFlushCache( hDstDS );
1341    if( CPLGetLastErrorType() != CE_None )
1342        bHasGotErr = TRUE;
1343    GDALClose( hDstDS );
1344
1345    setMapInMaps(outputs,"Result","value",pszDstFilename);
1346   
1347    CPLFree( pszDstFilename );
1348    CSLDestroy( papszSrcFiles );
1349    CSLDestroy( papszWarpOptions );
1350    CSLDestroy( papszTO );
1351
1352    GDALDumpOpenDatasets( stderr );
1353
1354    GDALDestroyDriverManager();
1355   
1356#ifdef OGR_ENABLED
1357    if( hCutline != NULL )
1358        OGR_G_DestroyGeometry( (OGRGeometryH) hCutline );
1359    OGRCleanupAll();
1360#endif
1361
1362
1363/* -------------------------------------------------------------------- */
1364/*      Clean allocated arguments.                                      */
1365/* -------------------------------------------------------------------- */
1366    free(dataPath);
1367    free(tempPath);
1368
1369    tmpMap=getMapFromMaps(inputs,"te","value");
1370    if(tmpMap!=NULL){
1371      free(pszCutlineDSName);
1372    }
1373
1374    tmpMap=getMapFromMaps(inputs,"cwhere","value");
1375    if(tmpMap!=NULL){
1376      free(pszCWHERE);
1377    }
1378
1379    tmpMap=getMapFromMaps(inputs,"cl","value");
1380    if(tmpMap!=NULL){
1381      free(pszCLayer);
1382    }
1383
1384    tmpMap=getMapFromMaps(inputs,"csql","value");
1385    if(tmpMap!=NULL){
1386      free(pszCSQL);
1387    }
1388
1389    return (bHasGotErr) ? SERVICE_FAILED : SERVICE_SUCCEEDED;
1390}
1391
1392/************************************************************************/
1393/*                        GDALWarpCreateOutput()                        */
1394/*                                                                      */
1395/*      Create the output file based on various commandline options,    */
1396/*      and the input file.                                             */
1397/*      If there's just one source file, then *phTransformArg and       */
1398/*      *phSrcDS will be set, in order them to be reused by main        */
1399/*      function. This saves dataset re-opening, and above all transform*/
1400/*      recomputation, which can be expensive in the -tps case          */
1401/************************************************************************/
1402
1403static GDALDatasetH
1404GDALWarpCreateOutput( maps*& conf,char **papszSrcFiles, const char *pszFilename, 
1405                      const char *pszFormat, char **papszTO, 
1406                      char ***ppapszCreateOptions, GDALDataType eDT,
1407                      void ** phTransformArg,
1408                      GDALDatasetH* phSrcDS)
1409
1410
1411{
1412    GDALDriverH hDriver;
1413    GDALDatasetH hDstDS;
1414    void *hTransformArg;
1415    GDALColorTableH hCT = NULL;
1416    double dfWrkMinX=0, dfWrkMaxX=0, dfWrkMinY=0, dfWrkMaxY=0;
1417    double dfWrkResX=0, dfWrkResY=0;
1418    int nDstBandCount = 0;
1419    std::vector<GDALColorInterp> apeColorInterpretations;
1420
1421    *phTransformArg = NULL;
1422    *phSrcDS = NULL;
1423
1424/* -------------------------------------------------------------------- */
1425/*      Find the output driver.                                         */
1426/* -------------------------------------------------------------------- */
1427    hDriver = GDALGetDriverByName( pszFormat );
1428    if( hDriver == NULL 
1429        || GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL ) == NULL )
1430    {
1431        int     iDr;
1432        char tmp[2048];
1433        sprintf( tmp, "Output driver `%s' not recognised or does not support\n"
1434                 "direct output file creation.  The following format drivers are configured\n"
1435                 "and support direct output:\n", 
1436                 pszFormat );
1437
1438        for( iDr = 0; iDr < GDALGetDriverCount(); iDr++ )
1439        {
1440            GDALDriverH hDriver = GDALGetDriver(iDr);
1441
1442            if( GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL) != NULL )
1443            {
1444              char *_tmp=strdup(tmp);
1445              sprintf( tmp,"%s  %s: %s\n",
1446                       _tmp,
1447                       GDALGetDriverShortName( hDriver  ),
1448                       GDALGetDriverLongName( hDriver ) );
1449              free(_tmp);
1450            }
1451        }
1452        setMapInMaps(conf,"lenv","message",tmp);
1453        GDALExit( 1 );
1454    }
1455
1456/* -------------------------------------------------------------------- */
1457/*      For virtual output files, we have to set a special subclass     */
1458/*      of dataset to create.                                           */
1459/* -------------------------------------------------------------------- */
1460    if( bVRT )
1461        *ppapszCreateOptions = 
1462            CSLSetNameValue( *ppapszCreateOptions, "SUBCLASS", 
1463                             "VRTWarpedDataset" );
1464
1465/* -------------------------------------------------------------------- */
1466/*      Loop over all input files to collect extents.                   */
1467/* -------------------------------------------------------------------- */
1468    int     iSrc;
1469    char    *pszThisTargetSRS = (char*)CSLFetchNameValue( papszTO, "DST_SRS" );
1470    if( pszThisTargetSRS != NULL )
1471        pszThisTargetSRS = CPLStrdup( pszThisTargetSRS );
1472
1473    for( iSrc = 0; papszSrcFiles[iSrc] != NULL; iSrc++ )
1474    {
1475        GDALDatasetH hSrcDS;
1476        const char *pszThisSourceSRS = CSLFetchNameValue(papszTO,"SRC_SRS");
1477
1478        hSrcDS = GDALOpen( papszSrcFiles[iSrc], GA_ReadOnly );
1479        if( hSrcDS == NULL )
1480            GDALExit( 1 );
1481
1482/* -------------------------------------------------------------------- */
1483/*      Check that there's at least one raster band                     */
1484/* -------------------------------------------------------------------- */
1485        if ( GDALGetRasterCount(hSrcDS) == 0 )
1486        {
1487            fprintf(stderr, "Input file %s has no raster bands.\n", papszSrcFiles[iSrc] );
1488            GDALExit( 1 );
1489        }
1490
1491        if( eDT == GDT_Unknown )
1492            eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS,1));
1493
1494/* -------------------------------------------------------------------- */
1495/*      If we are processing the first file, and it has a color         */
1496/*      table, then we will copy it to the destination file.            */
1497/* -------------------------------------------------------------------- */
1498        if( iSrc == 0 )
1499        {
1500            nDstBandCount = GDALGetRasterCount(hSrcDS);
1501            hCT = GDALGetRasterColorTable( GDALGetRasterBand(hSrcDS,1) );
1502            if( hCT != NULL )
1503            {
1504                hCT = GDALCloneColorTable( hCT );
1505                if( !bQuiet )
1506                    printf( "Copying color table from %s to new file.\n", 
1507                            papszSrcFiles[iSrc] );
1508            }
1509
1510            for(int iBand = 0; iBand < nDstBandCount; iBand++)
1511            {
1512                apeColorInterpretations.push_back(
1513                    GDALGetRasterColorInterpretation(GDALGetRasterBand(hSrcDS,iBand+1)) );
1514            }
1515        }
1516
1517/* -------------------------------------------------------------------- */
1518/*      Get the sourcesrs from the dataset, if not set already.         */
1519/* -------------------------------------------------------------------- */
1520        if( pszThisSourceSRS == NULL )
1521        {
1522            const char *pszMethod = CSLFetchNameValue( papszTO, "METHOD" );
1523
1524            if( GDALGetProjectionRef( hSrcDS ) != NULL 
1525                && strlen(GDALGetProjectionRef( hSrcDS )) > 0
1526                && (pszMethod == NULL || EQUAL(pszMethod,"GEOTRANSFORM")) )
1527                pszThisSourceSRS = GDALGetProjectionRef( hSrcDS );
1528           
1529            else if( GDALGetGCPProjection( hSrcDS ) != NULL
1530                     && strlen(GDALGetGCPProjection(hSrcDS)) > 0 
1531                     && GDALGetGCPCount( hSrcDS ) > 1 
1532                     && (pszMethod == NULL || EQUALN(pszMethod,"GCP_",4)) )
1533                pszThisSourceSRS = GDALGetGCPProjection( hSrcDS );
1534            else if( pszMethod != NULL && EQUAL(pszMethod,"RPC") )
1535                pszThisSourceSRS = SRS_WKT_WGS84;
1536            else
1537                pszThisSourceSRS = "";
1538        }
1539
1540        if( pszThisTargetSRS == NULL )
1541            pszThisTargetSRS = CPLStrdup( pszThisSourceSRS );
1542       
1543/* -------------------------------------------------------------------- */
1544/*      Create a transformation object from the source to               */
1545/*      destination coordinate system.                                  */
1546/* -------------------------------------------------------------------- */
1547        hTransformArg = 
1548            GDALCreateGenImgProjTransformer2( hSrcDS, NULL, papszTO );
1549       
1550        if( hTransformArg == NULL )
1551        {
1552            CPLFree( pszThisTargetSRS );
1553            GDALClose( hSrcDS );
1554            return NULL;
1555        }
1556       
1557        GDALTransformerInfo* psInfo = (GDALTransformerInfo*)hTransformArg;
1558
1559/* -------------------------------------------------------------------- */
1560/*      Get approximate output definition.                              */
1561/* -------------------------------------------------------------------- */
1562        double adfThisGeoTransform[6];
1563        double adfExtent[4];
1564        int    nThisPixels, nThisLines;
1565
1566        if( GDALSuggestedWarpOutput2( hSrcDS, 
1567                                      psInfo->pfnTransform, hTransformArg, 
1568                                      adfThisGeoTransform, 
1569                                      &nThisPixels, &nThisLines, 
1570                                      adfExtent, 0 ) != CE_None )
1571        {
1572            CPLFree( pszThisTargetSRS );
1573            GDALClose( hSrcDS );
1574            return NULL;
1575        }
1576       
1577        if (CPLGetConfigOption( "CHECK_WITH_INVERT_PROJ", NULL ) == NULL)
1578        {
1579            double MinX = adfExtent[0];
1580            double MaxX = adfExtent[2];
1581            double MaxY = adfExtent[3];
1582            double MinY = adfExtent[1];
1583            int bSuccess = TRUE;
1584           
1585            /* Check that the the edges of the target image are in the validity area */
1586            /* of the target projection */
1587#define N_STEPS 20
1588            int i,j;
1589            for(i=0;i<=N_STEPS && bSuccess;i++)
1590            {
1591                for(j=0;j<=N_STEPS && bSuccess;j++)
1592                {
1593                    double dfRatioI = i * 1.0 / N_STEPS;
1594                    double dfRatioJ = j * 1.0 / N_STEPS;
1595                    double expected_x = (1 - dfRatioI) * MinX + dfRatioI * MaxX;
1596                    double expected_y = (1 - dfRatioJ) * MinY + dfRatioJ * MaxY;
1597                    double x = expected_x;
1598                    double y = expected_y;
1599                    double z = 0;
1600                    /* Target SRS coordinates to source image pixel coordinates */
1601                    if (!psInfo->pfnTransform(hTransformArg, TRUE, 1, &x, &y, &z, &bSuccess) || !bSuccess)
1602                        bSuccess = FALSE;
1603                    /* Source image pixel coordinates to target SRS coordinates */
1604                    if (!psInfo->pfnTransform(hTransformArg, FALSE, 1, &x, &y, &z, &bSuccess) || !bSuccess)
1605                        bSuccess = FALSE;
1606                    if (fabs(x - expected_x) > (MaxX - MinX) / nThisPixels ||
1607                        fabs(y - expected_y) > (MaxY - MinY) / nThisLines)
1608                        bSuccess = FALSE;
1609                }
1610            }
1611           
1612            /* If not, retry with CHECK_WITH_INVERT_PROJ=TRUE that forces ogrct.cpp */
1613            /* to check the consistency of each requested projection result with the */
1614            /* invert projection */
1615            if (!bSuccess)
1616            {
1617                CPLSetConfigOption( "CHECK_WITH_INVERT_PROJ", "TRUE" );
1618                CPLDebug("WARP", "Recompute out extent with CHECK_WITH_INVERT_PROJ=TRUE");
1619
1620                if( GDALSuggestedWarpOutput2( hSrcDS, 
1621                                      psInfo->pfnTransform, hTransformArg, 
1622                                      adfThisGeoTransform, 
1623                                      &nThisPixels, &nThisLines, 
1624                                      adfExtent, 0 ) != CE_None )
1625                {
1626                    CPLFree( pszThisTargetSRS );
1627                    GDALClose( hSrcDS );
1628                    return NULL;
1629                }
1630            }
1631        }
1632
1633/* -------------------------------------------------------------------- */
1634/*      Expand the working bounds to include this region, ensure the    */
1635/*      working resolution is no more than this resolution.             */
1636/* -------------------------------------------------------------------- */
1637        if( dfWrkMaxX == 0.0 && dfWrkMinX == 0.0 )
1638        {
1639            dfWrkMinX = adfExtent[0];
1640            dfWrkMaxX = adfExtent[2];
1641            dfWrkMaxY = adfExtent[3];
1642            dfWrkMinY = adfExtent[1];
1643            dfWrkResX = adfThisGeoTransform[1];
1644            dfWrkResY = ABS(adfThisGeoTransform[5]);
1645        }
1646        else
1647        {
1648            dfWrkMinX = MIN(dfWrkMinX,adfExtent[0]);
1649            dfWrkMaxX = MAX(dfWrkMaxX,adfExtent[2]);
1650            dfWrkMaxY = MAX(dfWrkMaxY,adfExtent[3]);
1651            dfWrkMinY = MIN(dfWrkMinY,adfExtent[1]);
1652            dfWrkResX = MIN(dfWrkResX,adfThisGeoTransform[1]);
1653            dfWrkResY = MIN(dfWrkResY,ABS(adfThisGeoTransform[5]));
1654        }
1655
1656        if (iSrc == 0 && papszSrcFiles[1] == NULL)
1657        {
1658            *phTransformArg = hTransformArg;
1659            *phSrcDS = hSrcDS;
1660        }
1661        else
1662        {
1663            GDALDestroyGenImgProjTransformer( hTransformArg );
1664            GDALClose( hSrcDS );
1665        }
1666    }
1667
1668/* -------------------------------------------------------------------- */
1669/*      Did we have any usable sources?                                 */
1670/* -------------------------------------------------------------------- */
1671    if( nDstBandCount == 0 )
1672    {
1673        CPLError( CE_Failure, CPLE_AppDefined,
1674                  "No usable source images." );
1675        CPLFree( pszThisTargetSRS );
1676        return NULL;
1677    }
1678
1679/* -------------------------------------------------------------------- */
1680/*      Turn the suggested region into a geotransform and suggested     */
1681/*      number of pixels and lines.                                     */
1682/* -------------------------------------------------------------------- */
1683    double adfDstGeoTransform[6];
1684    int nPixels, nLines;
1685
1686    adfDstGeoTransform[0] = dfWrkMinX;
1687    adfDstGeoTransform[1] = dfWrkResX;
1688    adfDstGeoTransform[2] = 0.0;
1689    adfDstGeoTransform[3] = dfWrkMaxY;
1690    adfDstGeoTransform[4] = 0.0;
1691    adfDstGeoTransform[5] = -1 * dfWrkResY;
1692
1693    nPixels = (int) ((dfWrkMaxX - dfWrkMinX) / dfWrkResX + 0.5);
1694    nLines = (int) ((dfWrkMaxY - dfWrkMinY) / dfWrkResY + 0.5);
1695
1696/* -------------------------------------------------------------------- */
1697/*      Did the user override some parameters?                          */
1698/* -------------------------------------------------------------------- */
1699    if( dfXRes != 0.0 && dfYRes != 0.0 )
1700    {
1701        if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 )
1702        {
1703            dfMinX = adfDstGeoTransform[0];
1704            dfMaxX = adfDstGeoTransform[0] + adfDstGeoTransform[1] * nPixels;
1705            dfMaxY = adfDstGeoTransform[3];
1706            dfMinY = adfDstGeoTransform[3] + adfDstGeoTransform[5] * nLines;
1707        }
1708       
1709        if ( bTargetAlignedPixels )
1710        {
1711            dfMinX = floor(dfMinX / dfXRes) * dfXRes;
1712            dfMaxX = ceil(dfMaxX / dfXRes) * dfXRes;
1713            dfMinY = floor(dfMinY / dfYRes) * dfYRes;
1714            dfMaxY = ceil(dfMaxY / dfYRes) * dfYRes;
1715        }
1716
1717        nPixels = (int) ((dfMaxX - dfMinX + (dfXRes/2.0)) / dfXRes);
1718        nLines = (int) ((dfMaxY - dfMinY + (dfYRes/2.0)) / dfYRes);
1719        adfDstGeoTransform[0] = dfMinX;
1720        adfDstGeoTransform[3] = dfMaxY;
1721        adfDstGeoTransform[1] = dfXRes;
1722        adfDstGeoTransform[5] = -dfYRes;
1723    }
1724
1725    else if( nForcePixels != 0 && nForceLines != 0 )
1726    {
1727        if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 )
1728        {
1729            dfMinX = dfWrkMinX;
1730            dfMaxX = dfWrkMaxX;
1731            dfMaxY = dfWrkMaxY;
1732            dfMinY = dfWrkMinY;
1733        }
1734
1735        dfXRes = (dfMaxX - dfMinX) / nForcePixels;
1736        dfYRes = (dfMaxY - dfMinY) / nForceLines;
1737
1738        adfDstGeoTransform[0] = dfMinX;
1739        adfDstGeoTransform[3] = dfMaxY;
1740        adfDstGeoTransform[1] = dfXRes;
1741        adfDstGeoTransform[5] = -dfYRes;
1742
1743        nPixels = nForcePixels;
1744        nLines = nForceLines;
1745    }
1746
1747    else if( nForcePixels != 0 )
1748    {
1749        if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 )
1750        {
1751            dfMinX = dfWrkMinX;
1752            dfMaxX = dfWrkMaxX;
1753            dfMaxY = dfWrkMaxY;
1754            dfMinY = dfWrkMinY;
1755        }
1756
1757        dfXRes = (dfMaxX - dfMinX) / nForcePixels;
1758        dfYRes = dfXRes;
1759
1760        adfDstGeoTransform[0] = dfMinX;
1761        adfDstGeoTransform[3] = dfMaxY;
1762        adfDstGeoTransform[1] = dfXRes;
1763        adfDstGeoTransform[5] = -dfYRes;
1764
1765        nPixels = nForcePixels;
1766        nLines = (int) ((dfMaxY - dfMinY + (dfYRes/2.0)) / dfYRes);
1767    }
1768
1769    else if( nForceLines != 0 )
1770    {
1771        if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 )
1772        {
1773            dfMinX = dfWrkMinX;
1774            dfMaxX = dfWrkMaxX;
1775            dfMaxY = dfWrkMaxY;
1776            dfMinY = dfWrkMinY;
1777        }
1778
1779        dfYRes = (dfMaxY - dfMinY) / nForceLines;
1780        dfXRes = dfYRes;
1781
1782        adfDstGeoTransform[0] = dfMinX;
1783        adfDstGeoTransform[3] = dfMaxY;
1784        adfDstGeoTransform[1] = dfXRes;
1785        adfDstGeoTransform[5] = -dfYRes;
1786
1787        nPixels = (int) ((dfMaxX - dfMinX + (dfXRes/2.0)) / dfXRes);
1788        nLines = nForceLines;
1789    }
1790
1791    else if( dfMinX != 0.0 || dfMinY != 0.0 || dfMaxX != 0.0 || dfMaxY != 0.0 )
1792    {
1793        dfXRes = adfDstGeoTransform[1];
1794        dfYRes = fabs(adfDstGeoTransform[5]);
1795
1796        nPixels = (int) ((dfMaxX - dfMinX + (dfXRes/2.0)) / dfXRes);
1797        nLines = (int) ((dfMaxY - dfMinY + (dfYRes/2.0)) / dfYRes);
1798
1799        dfXRes = (dfMaxX - dfMinX) / nPixels;
1800        dfYRes = (dfMaxY - dfMinY) / nLines;
1801
1802        adfDstGeoTransform[0] = dfMinX;
1803        adfDstGeoTransform[3] = dfMaxY;
1804        adfDstGeoTransform[1] = dfXRes;
1805        adfDstGeoTransform[5] = -dfYRes;
1806    }
1807
1808/* -------------------------------------------------------------------- */
1809/*      Do we want to generate an alpha band in the output file?        */
1810/* -------------------------------------------------------------------- */
1811    if( bEnableSrcAlpha )
1812        nDstBandCount--;
1813
1814    if( bEnableDstAlpha )
1815        nDstBandCount++;
1816
1817/* -------------------------------------------------------------------- */
1818/*      Create the output file.                                         */
1819/* -------------------------------------------------------------------- */
1820    if( !bQuiet )
1821        printf( "Creating output file that is %dP x %dL.\n", nPixels, nLines );
1822
1823    hDstDS = GDALCreate( hDriver, pszFilename, nPixels, nLines, 
1824                         nDstBandCount, eDT, *ppapszCreateOptions );
1825   
1826    if( hDstDS == NULL )
1827    {
1828        CPLFree( pszThisTargetSRS );
1829        return NULL;
1830    }
1831
1832/* -------------------------------------------------------------------- */
1833/*      Write out the projection definition.                            */
1834/* -------------------------------------------------------------------- */
1835    GDALSetProjection( hDstDS, pszThisTargetSRS );
1836    GDALSetGeoTransform( hDstDS, adfDstGeoTransform );
1837
1838    if (*phTransformArg != NULL)
1839        GDALSetGenImgProjTransformerDstGeoTransform( *phTransformArg, adfDstGeoTransform);
1840
1841/* -------------------------------------------------------------------- */
1842/*      Try to set color interpretation of source bands to target       */
1843/*      dataset.                                                        */
1844/*      FIXME? We should likely do that for other drivers than VRT      */
1845/*      but it might create spurious .aux.xml files (at least with HFA, */
1846/*      and netCDF)                                                     */
1847/* -------------------------------------------------------------------- */
1848    if( bVRT )
1849    {
1850        int nBandsToCopy = (int)apeColorInterpretations.size();
1851        if ( bEnableSrcAlpha )
1852            nBandsToCopy --;
1853        for(int iBand = 0; iBand < nBandsToCopy; iBand++)
1854        {
1855            GDALSetRasterColorInterpretation(
1856                GDALGetRasterBand( hDstDS, iBand + 1 ),
1857                apeColorInterpretations[iBand] );
1858        }
1859    }
1860   
1861/* -------------------------------------------------------------------- */
1862/*      Try to set color interpretation of output file alpha band.      */
1863/* -------------------------------------------------------------------- */
1864    if( bEnableDstAlpha )
1865    {
1866        GDALSetRasterColorInterpretation( 
1867            GDALGetRasterBand( hDstDS, nDstBandCount ), 
1868            GCI_AlphaBand );
1869    }
1870
1871/* -------------------------------------------------------------------- */
1872/*      Copy the color table, if required.                              */
1873/* -------------------------------------------------------------------- */
1874    if( hCT != NULL )
1875    {
1876        GDALSetRasterColorTable( GDALGetRasterBand(hDstDS,1), hCT );
1877        GDALDestroyColorTable( hCT );
1878    }
1879
1880    CPLFree( pszThisTargetSRS );
1881    return hDstDS;
1882}
1883
1884/************************************************************************/
1885/*                      GeoTransform_Transformer()                      */
1886/*                                                                      */
1887/*      Convert points from georef coordinates to pixel/line based      */
1888/*      on a geotransform.                                              */
1889/************************************************************************/
1890
1891class CutlineTransformer : public OGRCoordinateTransformation
1892{
1893public:
1894
1895    void         *hSrcImageTransformer;
1896
1897    virtual OGRSpatialReference *GetSourceCS() { return NULL; }
1898    virtual OGRSpatialReference *GetTargetCS() { return NULL; }
1899
1900    virtual int Transform( int nCount, 
1901                           double *x, double *y, double *z = NULL ) {
1902        int nResult;
1903
1904        int *pabSuccess = (int *) CPLCalloc(sizeof(int),nCount);
1905        nResult = TransformEx( nCount, x, y, z, pabSuccess );
1906        CPLFree( pabSuccess );
1907
1908        return nResult;
1909    }
1910
1911    virtual int TransformEx( int nCount, 
1912                             double *x, double *y, double *z = NULL,
1913                             int *pabSuccess = NULL ) {
1914        return GDALGenImgProjTransform( hSrcImageTransformer, TRUE, 
1915                                        nCount, x, y, z, pabSuccess );
1916    }
1917};
1918
1919
1920/************************************************************************/
1921/*                            LoadCutline()                             */
1922/*                                                                      */
1923/*      Load blend cutline from OGR datasource.                         */
1924/************************************************************************/
1925
1926static void
1927LoadCutline( const char *pszCutlineDSName, const char *pszCLayer, 
1928             const char *pszCWHERE, const char *pszCSQL, 
1929             void **phCutlineRet )
1930
1931{
1932#ifndef OGR_ENABLED
1933    CPLError( CE_Failure, CPLE_AppDefined, 
1934              "Request to load a cutline failed, this build does not support OGR features.\n" );
1935    GDALExit( 1 );
1936#else // def OGR_ENABLED
1937    OGRRegisterAll();
1938
1939/* -------------------------------------------------------------------- */
1940/*      Open source vector dataset.                                     */
1941/* -------------------------------------------------------------------- */
1942    OGRDataSourceH hSrcDS;
1943
1944    hSrcDS = OGROpen( pszCutlineDSName, FALSE, NULL );
1945    if( hSrcDS == NULL )
1946        GDALExit( 1 );
1947
1948/* -------------------------------------------------------------------- */
1949/*      Get the source layer                                            */
1950/* -------------------------------------------------------------------- */
1951    OGRLayerH hLayer = NULL;
1952
1953    if( pszCSQL != NULL )
1954        hLayer = OGR_DS_ExecuteSQL( hSrcDS, pszCSQL, NULL, NULL ); 
1955    else if( pszCLayer != NULL )
1956        hLayer = OGR_DS_GetLayerByName( hSrcDS, pszCLayer );
1957    else
1958        hLayer = OGR_DS_GetLayer( hSrcDS, 0 );
1959
1960    if( hLayer == NULL )
1961    {
1962        fprintf( stderr, "Failed to identify source layer from datasource.\n" );
1963        GDALExit( 1 );
1964    }
1965
1966/* -------------------------------------------------------------------- */
1967/*      Apply WHERE clause if there is one.                             */
1968/* -------------------------------------------------------------------- */
1969    if( pszCWHERE != NULL )
1970        OGR_L_SetAttributeFilter( hLayer, pszCWHERE );
1971
1972/* -------------------------------------------------------------------- */
1973/*      Collect the geometries from this layer, and build list of       */
1974/*      burn values.                                                    */
1975/* -------------------------------------------------------------------- */
1976    OGRFeatureH hFeat;
1977    OGRGeometryH hMultiPolygon = OGR_G_CreateGeometry( wkbMultiPolygon );
1978
1979    OGR_L_ResetReading( hLayer );
1980   
1981    while( (hFeat = OGR_L_GetNextFeature( hLayer )) != NULL )
1982    {
1983        OGRGeometryH hGeom = OGR_F_GetGeometryRef(hFeat);
1984
1985        if( hGeom == NULL )
1986        {
1987            fprintf( stderr, "ERROR: Cutline feature without a geometry.\n" );
1988            GDALExit( 1 );
1989        }
1990       
1991        OGRwkbGeometryType eType = wkbFlatten(OGR_G_GetGeometryType( hGeom ));
1992
1993        if( eType == wkbPolygon )
1994            OGR_G_AddGeometry( hMultiPolygon, hGeom );
1995        else if( eType == wkbMultiPolygon )
1996        {
1997            int iGeom;
1998
1999            for( iGeom = 0; iGeom < OGR_G_GetGeometryCount( hGeom ); iGeom++ )
2000            {
2001                OGR_G_AddGeometry( hMultiPolygon, 
2002                                   OGR_G_GetGeometryRef(hGeom,iGeom) );
2003            }
2004        }
2005        else
2006        {
2007            fprintf( stderr, "ERROR: Cutline not of polygon type.\n" );
2008            GDALExit( 1 );
2009        }
2010
2011        OGR_F_Destroy( hFeat );
2012    }
2013
2014    if( OGR_G_GetGeometryCount( hMultiPolygon ) == 0 )
2015    {
2016        fprintf( stderr, "ERROR: Did not get any cutline features.\n" );
2017        GDALExit( 1 );
2018    }
2019
2020/* -------------------------------------------------------------------- */
2021/*      Ensure the coordinate system gets set on the geometry.          */
2022/* -------------------------------------------------------------------- */
2023    OGR_G_AssignSpatialReference(
2024        hMultiPolygon, OGR_L_GetSpatialRef(hLayer) );
2025
2026    *phCutlineRet = (void *) hMultiPolygon;
2027
2028/* -------------------------------------------------------------------- */
2029/*      Cleanup                                                         */
2030/* -------------------------------------------------------------------- */
2031    if( pszCSQL != NULL )
2032        OGR_DS_ReleaseResultSet( hSrcDS, hLayer );
2033
2034    OGR_DS_Destroy( hSrcDS );
2035#endif
2036}
2037
2038/************************************************************************/
2039/*                      TransformCutlineToSource()                      */
2040/*                                                                      */
2041/*      Transform cutline from its SRS to source pixel/line coordinates.*/
2042/************************************************************************/
2043static void
2044TransformCutlineToSource( GDALDatasetH hSrcDS, void *hCutline,
2045                          char ***ppapszWarpOptions, char **papszTO_In )
2046
2047{
2048#ifdef OGR_ENABLED
2049    OGRGeometryH hMultiPolygon = OGR_G_Clone( (OGRGeometryH) hCutline );
2050    char **papszTO = CSLDuplicate( papszTO_In );
2051
2052/* -------------------------------------------------------------------- */
2053/*      Checkout that SRS are the same.                                 */
2054/* -------------------------------------------------------------------- */
2055    OGRSpatialReferenceH  hRasterSRS = NULL;
2056    const char *pszProjection = NULL;
2057
2058    if( GDALGetProjectionRef( hSrcDS ) != NULL 
2059        && strlen(GDALGetProjectionRef( hSrcDS )) > 0 )
2060        pszProjection = GDALGetProjectionRef( hSrcDS );
2061    else if( GDALGetGCPProjection( hSrcDS ) != NULL )
2062        pszProjection = GDALGetGCPProjection( hSrcDS );
2063
2064    if( pszProjection != NULL )
2065    {
2066        hRasterSRS = OSRNewSpatialReference(NULL);
2067        if( OSRImportFromWkt( hRasterSRS, (char **)&pszProjection ) != CE_None )
2068        {
2069            OSRDestroySpatialReference(hRasterSRS);
2070            hRasterSRS = NULL;
2071        }
2072    }
2073
2074    OGRSpatialReferenceH hCutlineSRS = OGR_G_GetSpatialReference( hMultiPolygon );
2075    if( hRasterSRS != NULL && hCutlineSRS != NULL )
2076    {
2077        /* ok, we will reproject */
2078    }
2079    else if( hRasterSRS != NULL && hCutlineSRS == NULL )
2080    {
2081        fprintf(stderr,
2082                "Warning : the source raster dataset has a SRS, but the cutline features\n"
2083                "not.  We assume that the cutline coordinates are expressed in the destination SRS.\n"
2084                "If not, cutline results may be incorrect.\n");
2085    }
2086    else if( hRasterSRS == NULL && hCutlineSRS != NULL )
2087    {
2088        fprintf(stderr,
2089                "Warning : the input vector layer has a SRS, but the source raster dataset does not.\n"
2090                "Cutline results may be incorrect.\n");
2091    }
2092
2093    if( hRasterSRS != NULL )
2094        OSRDestroySpatialReference(hRasterSRS);
2095
2096/* -------------------------------------------------------------------- */
2097/*      Extract the cutline SRS WKT.                                    */
2098/* -------------------------------------------------------------------- */
2099    if( hCutlineSRS != NULL )
2100    {
2101        char *pszCutlineSRS_WKT = NULL;
2102
2103        OSRExportToWkt( hCutlineSRS, &pszCutlineSRS_WKT );
2104        papszTO = CSLSetNameValue( papszTO, "DST_SRS", pszCutlineSRS_WKT );
2105        CPLFree( pszCutlineSRS_WKT );
2106    }
2107
2108/* -------------------------------------------------------------------- */
2109/*      It may be unwise to let the mask geometry be re-wrapped by      */
2110/*      the CENTER_LONG machinery as this can easily screw up world     */
2111/*      spanning masks and invert the mask topology.                    */
2112/* -------------------------------------------------------------------- */
2113    papszTO = CSLSetNameValue( papszTO, "INSERT_CENTER_LONG", "FALSE" );
2114
2115/* -------------------------------------------------------------------- */
2116/*      Transform the geometry to pixel/line coordinates.               */
2117/* -------------------------------------------------------------------- */
2118    CutlineTransformer oTransformer;
2119
2120    /* The cutline transformer will *invert* the hSrcImageTransformer */
2121    /* so it will convert from the cutline SRS to the source pixel/line */
2122    /* coordinates */
2123    oTransformer.hSrcImageTransformer = 
2124        GDALCreateGenImgProjTransformer2( hSrcDS, NULL, papszTO );
2125
2126    CSLDestroy( papszTO );
2127
2128    if( oTransformer.hSrcImageTransformer == NULL )
2129        GDALExit( 1 );
2130
2131    OGR_G_Transform( hMultiPolygon, 
2132                     (OGRCoordinateTransformationH) &oTransformer );
2133
2134    GDALDestroyGenImgProjTransformer( oTransformer.hSrcImageTransformer );
2135
2136/* -------------------------------------------------------------------- */
2137/*      Convert aggregate geometry into WKT.                            */
2138/* -------------------------------------------------------------------- */
2139    char *pszWKT = NULL;
2140
2141    OGR_G_ExportToWkt( hMultiPolygon, &pszWKT );
2142    OGR_G_DestroyGeometry( hMultiPolygon );
2143
2144    *ppapszWarpOptions = CSLSetNameValue( *ppapszWarpOptions, 
2145                                          "CUTLINE", pszWKT );
2146    CPLFree( pszWKT );
2147#endif
2148}
2149
2150}
Note: See TracBrowser for help on using the repository browser.

Search

Context Navigation

ZOO Sponsors

http://www.zoo-project.org/trac/chrome/site/img/geolabs-logo.pnghttp://www.zoo-project.org/trac/chrome/site/img/neogeo-logo.png http://www.zoo-project.org/trac/chrome/site/img/apptech-logo.png http://www.zoo-project.org/trac/chrome/site/img/3liz-logo.png http://www.zoo-project.org/trac/chrome/site/img/gateway-logo.png

Become a sponsor !

Knowledge partners

http://www.zoo-project.org/trac/chrome/site/img/ocu-logo.png http://www.zoo-project.org/trac/chrome/site/img/gucas-logo.png http://www.zoo-project.org/trac/chrome/site/img/polimi-logo.png http://www.zoo-project.org/trac/chrome/site/img/fem-logo.png http://www.zoo-project.org/trac/chrome/site/img/supsi-logo.png http://www.zoo-project.org/trac/chrome/site/img/cumtb-logo.png

Become a knowledge partner

Related links

http://zoo-project.org/img/ogclogo.png http://zoo-project.org/img/osgeologo.png