source: trunk/zoo-project/zoo-services/gdal/translate/service.c @ 369

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

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

  • Property svn:eol-style set to native
  • Property svn:mime-type set to text/x-csrc
File size: 29.7 KB
Line 
1/******************************************************************************
2 * $Id$
3 *
4 * Project:  GDAL Utilities
5 * Purpose:  GDAL Image Translator Program
6 * Author:   Frank Warmerdam, warmerdam@pobox.com
7 *
8 ******************************************************************************
9 * Copyright (c) 1998, 2002, Frank Warmerdam
10 *
11 * Permission is hereby granted, free of charge, to any person obtaining a
12 * copy of this software and associated documentation files (the "Software"),
13 * to deal in the Software without restriction, including without limitation
14 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 * and/or sell copies of the Software, and to permit persons to whom the
16 * Software is furnished to do so, subject to the following conditions:
17 *
18 * The above copyright notice and this permission notice shall be included
19 * in all copies or substantial portions of the Software.
20 *
21 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 * DEALINGS IN THE SOFTWARE.
28 ****************************************************************************/
29
30#include "cpl_vsi.h"
31#include "cpl_conv.h"
32#include "cpl_string.h"
33#include "gdal_priv.h"
34#include "ogr_spatialref.h"
35#include "vrt/vrtdataset.h"
36
37#include "service.h"
38
39CPL_CVSID("$Id$");
40
41extern "C" {
42
43  static void AttachMetadata( GDALDatasetH, char ** );
44  static int bSubCall = FALSE;
45
46  /************************************************************************/
47  /*                          Gdal_Translate()                            */
48  /************************************************************************/
49
50#ifdef WIN32
51__declspec(dllexport)
52#endif
53  int Gdal_Translate(maps*& conf,maps*& inputs,maps*& outputs)
54  {
55
56
57    GDALDatasetH        hDataset, hOutDS;
58    int                 i;
59    int                 nRasterXSize, nRasterYSize;
60    const char          *pszSource=NULL, *pszDest=NULL, *pszFormat = "GTiff";
61    GDALDriverH         hDriver;
62    int                 *panBandList = NULL, nBandCount = 0, bDefBands = TRUE;
63    double              adfGeoTransform[6];
64    GDALDataType        eOutputType = GDT_Unknown;
65    int                 nOXSize = 0, nOYSize = 0;
66    char                *pszOXSize=NULL, *pszOYSize=NULL;
67    char                **papszCreateOptions = NULL;
68    int                 anSrcWin[4], bStrict = FALSE;
69    const char          *pszProjection;
70    int                 bScale = FALSE, bHaveScaleSrc = FALSE;
71    double              dfScaleSrcMin=0.0, dfScaleSrcMax=255.0;
72    double              dfScaleDstMin=0.0, dfScaleDstMax=255.0;
73    double              dfULX, dfULY, dfLRX, dfLRY;
74    char                **papszMetadataOptions = NULL;
75    char                *pszOutputSRS = NULL;
76    int                 bQuiet = TRUE, bGotBounds = FALSE;
77    GDALProgressFunc    pfnProgress = GDALDummyProgress;
78    int                 nGCPCount = 0;
79    GDAL_GCP            *pasGCPs = NULL;
80    int                 iSrcFileArg = -1, iDstFileArg = -1;
81    int                 bCopySubDatasets = FALSE;
82    double              adfULLR[4] = { 0,0,0,0 };
83    int                 bSetNoData = FALSE;
84    double              dfNoDataReal = 0.0;
85    int                 nRGBExpand = 0;
86
87    anSrcWin[0] = 0;
88    anSrcWin[1] = 0;
89    anSrcWin[2] = 0;
90    anSrcWin[3] = 0;
91
92    dfULX = dfULY = dfLRX = dfLRY = 0.0;
93
94    /* ----------------------------------------------------------------- */
95    /*      Register standard GDAL drivers, and process generic GDAL     */
96    /* ----------------------------------------------------------------- */
97    GDALAllRegister();
98    /* ----------------------------------------------------------------- */
99    /* Extract Format, InputDSN, OutputDSN parameters                    */
100    /* ----------------------------------------------------------------- */
101
102    map* tmpMap=NULL;
103
104    char dataPath[1024];
105    tmpMap=getMapFromMaps(conf,"main","dataPath");
106    if(tmpMap!=NULL)
107      sprintf(dataPath,"%s",tmpMap->value);
108    tmpMap=NULL;
109
110    char tempPath[1024];
111    tmpMap=getMapFromMaps(conf,"main","tmpPath");
112    if(tmpMap!=NULL){
113      sprintf(tempPath,"%s",tmpMap->value);
114    }
115    tmpMap=NULL;
116
117    tmpMap=getMapFromMaps(inputs,"Format","value");
118    if(tmpMap!=NULL){
119      pszFormat=tmpMap->value;
120    }
121    tmpMap=NULL;
122    tmpMap=getMapFromMaps(inputs,"InputDSN","value");
123    if(tmpMap!=NULL){
124      pszSource=(char*)CPLMalloc(sizeof(char)*(strlen(dataPath)+strlen(tmpMap->value)+4));
125      sprintf((char*)pszSource,"%s/%s",dataPath,tmpMap->value);
126    }
127    tmpMap=NULL;
128    tmpMap=getMapFromMaps(inputs,"OutputDSN","value");
129    if(tmpMap!=NULL){
130      pszDest=(char*)CPLMalloc(sizeof(char)*(strlen(tempPath)+strlen(tmpMap->value)+4));
131      char *ext=new char[4];
132      ext="tif";
133      if(strncasecmp(pszFormat,"AAIGRID",7)==0)
134        ext="csv";
135      else 
136        if(strncasecmp(pszFormat,"PNG",3)==0)
137          ext="png";
138        else
139          if(strncasecmp(pszFormat,"GIF",3)==0)
140            ext="gif";
141          else
142            if(strncasecmp(pszFormat,"JPEG",4)==0)
143              ext="jpg";
144      sprintf((char*)pszDest,"%s/%s.%s",tempPath,tmpMap->value,ext);
145      fprintf(stderr,"DEBUG pszDest : %s\n",pszDest);
146    }
147    tmpMap=NULL;
148    tmpMap=getMapFromMaps(inputs,"ProjWin","value");
149    if(tmpMap!=NULL){
150      char *tmp=tmpMap->value;
151      char *t=strtok(tmp,",");
152      int cnt=0;
153      while(t!=NULL){
154        switch(cnt){
155        case 0:
156          dfULX = atof(t);
157          break;
158        case 1:
159          dfULY = atof(t);
160          break;
161        case 2:
162          dfLRX = atof(t);
163          break;
164        case 3:
165          dfLRY = atof(t);
166          break;
167        }
168        t=strtok(NULL,",");
169        cnt++;
170      }
171    }
172
173    tmpMap=NULL;
174    tmpMap=getMapFromMaps(inputs,"GCP","value");
175    if(tmpMap!=NULL){
176      map* length=getMapFromMaps(inputs,"GCP","length");
177      int len=0;
178      if(length){
179        len=atoi(length->value);
180        int i;
181        maps* currentMaps=getMaps(inputs,"GCP");
182        for(i=0;i<len;i++){
183          char* endptr = NULL;
184          /* -gcp pixel line easting northing [elev] */
185         
186          nGCPCount++;
187          pasGCPs = (GDAL_GCP *) 
188            realloc( pasGCPs, sizeof(GDAL_GCP) * nGCPCount );
189          GDALInitGCPs( 1, pasGCPs + nGCPCount - 1 );
190
191          map* currentMap=getMapArray(currentMaps->content,"value",i);
192
193          char* tmpV=strdup(currentMap->value);
194          char *res=strtok(tmpV,",");
195          int j=0;
196          while(res!=NULL){
197            switch(j){
198            case 0:
199              pasGCPs[nGCPCount-1].dfGCPPixel = CPLAtofM(res);
200              break;
201            case 1:
202              pasGCPs[nGCPCount-1].dfGCPLine = CPLAtofM(res);
203              break;
204            case 2:
205              pasGCPs[nGCPCount-1].dfGCPX = CPLAtofM(res);
206              break;
207            case 3:
208              pasGCPs[nGCPCount-1].dfGCPY = CPLAtofM(res);
209              break;
210            case 4:
211              if(res!=NULL && (strtod(res, &endptr) != 0.0 || res[0] == '0'))
212                if (endptr && *endptr == 0)
213                  pasGCPs[nGCPCount-1].dfGCPZ = CPLAtofM(res);
214              break;
215            }
216            res=strtok(NULL,",");
217            j++;
218          }
219        }
220      }else{
221        char* endptr = NULL;
222        /* -gcp pixel line easting northing [elev] */
223       
224        nGCPCount++;
225        pasGCPs = (GDAL_GCP *) 
226          realloc( pasGCPs, sizeof(GDAL_GCP) * nGCPCount );
227        GDALInitGCPs( 1, pasGCPs + nGCPCount - 1 );
228       
229        char* tmpV=strdup(tmpMap->value);
230        char *res=strtok(tmpV,",");
231        int j=0;
232        while(res!=NULL){
233          switch(j){
234          case 0:
235            pasGCPs[nGCPCount-1].dfGCPPixel = CPLAtofM(res);
236            break;
237          case 1:
238            pasGCPs[nGCPCount-1].dfGCPLine = CPLAtofM(res);
239            break;
240          case 2:
241            pasGCPs[nGCPCount-1].dfGCPX = CPLAtofM(res);
242            break;
243          case 3:
244            pasGCPs[nGCPCount-1].dfGCPY = CPLAtofM(res);
245            break;
246          case 4:
247            if(res!=NULL && (CPLStrtod(res, &endptr) != 0.0 || res[0] == '0'))
248              if (endptr && *endptr == 0)
249                pasGCPs[nGCPCount-1].dfGCPZ = CPLAtofM(res);
250            break;
251          }         
252          res=strtok(NULL,",");
253          j++;
254        }
255      }
256    }
257
258    tmpMap=NULL;
259    tmpMap=getMapFromMaps(inputs,"SRS","value");
260    if(tmpMap!=NULL){
261      OGRSpatialReference oOutputSRS;
262      if( oOutputSRS.SetFromUserInput( tmpMap->value ) != OGRERR_NONE )
263        {
264          char *msg=(char*)CPLMalloc(100*sizeof(char));
265          sprintf( msg, "Failed to process SRS definition: %s\n", 
266                   tmpMap->value );
267          setMapInMaps(conf,"lenv","message",msg);
268          /**
269           * Avoiding GDALDestroyDriverManager() call
270           */
271          return SERVICE_FAILED;
272        }
273      oOutputSRS.exportToWkt( &pszOutputSRS );
274    }
275    tmpMap=NULL;
276    tmpMap=getMapFromMaps(inputs,"Type","value");
277    if(tmpMap!=NULL){
278      int       iType;
279     
280      for( iType = 1; iType < GDT_TypeCount; iType++ )
281        {
282          if( GDALGetDataTypeName((GDALDataType)iType) != NULL
283              && EQUAL(GDALGetDataTypeName((GDALDataType)iType),
284                       tmpMap->value) )
285            {
286              eOutputType = (GDALDataType) iType;
287            }
288        }
289     
290      if( eOutputType == GDT_Unknown )
291        {
292          printf( "Unknown output pixel type: %s\n", tmpMap->value );
293          /**
294           * Avoiding GDALDestroyDriverManager() call
295           */
296          exit( 2 );
297        }
298    }
299
300    if( pszDest == NULL ){
301        fprintf(stderr,"exit line 416");
302        fflush(stderr);
303        /**
304         * Avoiding GDALDestroyDriverManager() call
305         */
306        exit( 10 );
307      }
308
309    if ( strcmp(pszSource, pszDest) == 0)
310      {
311        fprintf(stderr, "Source and destination datasets must be different.\n");
312        fflush(stderr);
313        /**
314         * Avoiding GDALDestroyDriverManager() call
315         */
316        exit( 1 );
317      }
318
319    /* ----------------------------------------------------------------- */
320    /*      Attempt to open source file.                                 */
321    /* ----------------------------------------------------------------- */
322
323    hDataset = GDALOpenShared( pszSource, GA_ReadOnly );
324   
325    if( hDataset == NULL ){
326      char *msg=(char*) CPLMalloc(1024*sizeof(char));
327      sprintf( msg,
328               "GDALOpen failed - %d\n%s\n",
329               CPLGetLastErrorNo(), CPLGetLastErrorMsg() );
330      setMapInMaps(conf,"lenv","message",msg);
331      return SERVICE_FAILED;
332    }
333
334    /* ----------------------------------------------------------------- */
335    /*      Handle subdatasets.                                          */
336    /* ----------------------------------------------------------------- */
337    if( !bCopySubDatasets
338        && CSLCount(GDALGetMetadata( hDataset, "SUBDATASETS" )) > 0 
339        && GDALGetRasterCount(hDataset) == 0 )
340      {
341        char *msg=(char*) CPLMalloc(1024*sizeof(char));
342        sprintf( msg,
343                 "Input file contains subdatasets. Please, select one of them for reading.\n" );
344        setMapInMaps(conf,"lenv","message",msg);
345        GDALClose( hDataset );
346        /**
347         * Avoiding GDALDestroyDriverManager() call
348         */
349        exit( 1 );
350      }
351
352    if( CSLCount(GDALGetMetadata( hDataset, "SUBDATASETS" )) > 0 
353        && bCopySubDatasets )
354      {
355        char **papszSubdatasets = GDALGetMetadata(hDataset,"SUBDATASETS");
356        char *pszSubDest = (char *) CPLMalloc(strlen(pszDest)+32);
357        int i;
358        int bOldSubCall = bSubCall;
359       
360        //argv[iDstFileArg] = pszSubDest;
361        bSubCall = TRUE;
362        for( i = 0; papszSubdatasets[i] != NULL; i += 2 )
363          {
364            //argv[iSrcFileArg] = strstr(papszSubdatasets[i],"=")+1;
365            sprintf( pszSubDest, "%s%d", pszDest, i/2 + 1 );
366            /*if( ProxyMain( argc, argv ) != 0 )
367              break;*/
368          }
369       
370        bSubCall = bOldSubCall;
371        free( pszSubDest );
372
373        GDALClose( hDataset );
374
375        if( !bSubCall )
376          {
377            GDALDumpOpenDatasets( stderr );
378            fflush(stderr);
379            /**
380             * Avoiding GDALDestroyDriverManager() call
381             */
382          }
383        return 1;
384      }
385
386    /* ----------------------------------------------------------------- */
387    /*      Collect some information from the source file.               */
388    /* ----------------------------------------------------------------- */
389    nRasterXSize = GDALGetRasterXSize( hDataset );
390    nRasterYSize = GDALGetRasterYSize( hDataset );
391
392    if( !bQuiet )
393      fprintf( stderr, "Input file size is %d, %d\n", nRasterXSize, nRasterYSize );
394
395    if( anSrcWin[2] == 0 && anSrcWin[3] == 0 ){
396        anSrcWin[2] = nRasterXSize;
397        anSrcWin[3] = nRasterYSize;
398      }
399
400    /* ----------------------------------------------------------------- */
401    /*  Build band list to translate                                     */
402    /* ----------------------------------------------------------------- */
403    if( nBandCount == 0 ){
404        nBandCount = GDALGetRasterCount( hDataset );
405        if( nBandCount == 0 ){
406            fprintf( stderr, "Input file has no bands, and so cannot be translated.\n" );
407            fflush(stderr);
408            /**
409             * Avoiding GDALDestroyDriverManager() call
410             */
411            exit(1 );
412          }
413
414        panBandList = (int *) CPLMalloc(sizeof(int)*nBandCount);
415        for( i = 0; i < nBandCount; i++ )
416          panBandList[i] = i+1;
417      }
418    else
419      {
420        for( i = 0; i < nBandCount; i++ )
421          {
422            if( panBandList[i] < 1 || panBandList[i] > GDALGetRasterCount(hDataset) )
423              {
424                fprintf( stderr, 
425                         "Band %d requested, but only bands 1 to %d available.\n",
426                         panBandList[i], GDALGetRasterCount(hDataset) );
427                fflush(stderr);
428                /**
429                 * Avoiding GDALDestroyDriverManager() call
430                 */
431                exit( 2 );
432              }
433          }
434
435        if( nBandCount != GDALGetRasterCount( hDataset ) )
436          bDefBands = FALSE;
437      }
438
439    /* ----------------------------------------------------------------- */
440    /*   Compute the source window from the projected source window      */
441    /*   if the projected coordinates were provided.  Note that the      */
442    /*   projected coordinates are in ulx, uly, lrx, lry format,         */
443    /*   while the anSrcWin is xoff, yoff, xsize, ysize with the         */
444    /*   xoff,yoff being the ulx, uly in pixel/line.                     */
445    /* ----------------------------------------------------------------- */
446    if( dfULX != 0.0 || dfULY != 0.0 
447        || dfLRX != 0.0 || dfLRY != 0.0 )
448      {
449        double  adfGeoTransform[6];
450
451        GDALGetGeoTransform( hDataset, adfGeoTransform );
452
453        if( adfGeoTransform[2] != 0.0 || adfGeoTransform[4] != 0.0 ){
454            fprintf( stderr, 
455                     "The -projwin option was used, but the geotransform is\n"
456                     "rotated.  This configuration is not supported.\n" );
457            GDALClose( hDataset );
458            free( panBandList );
459            fflush(stderr);
460            /**
461             * Avoiding GDALDestroyDriverManager() call
462             */
463            exit( 1 );
464          }
465
466        anSrcWin[0] = (int) 
467          ((dfULX - adfGeoTransform[0]) / adfGeoTransform[1] + 0.001);
468        anSrcWin[1] = (int) 
469          ((dfULY - adfGeoTransform[3]) / adfGeoTransform[5] + 0.001);
470
471        anSrcWin[2] = (int) ((dfLRX - dfULX) / adfGeoTransform[1] + 0.5);
472        anSrcWin[3] = (int) ((dfLRY - dfULY) / adfGeoTransform[5] + 0.5);
473
474        if( !bQuiet )
475          fprintf( stdout, 
476                   "Computed -srcwin %d %d %d %d from projected window.\n",
477                   anSrcWin[0], 
478                   anSrcWin[1], 
479                   anSrcWin[2], 
480                   anSrcWin[3] );
481       
482        if( anSrcWin[0] < 0 || anSrcWin[1] < 0 
483            || anSrcWin[0] + anSrcWin[2] > GDALGetRasterXSize(hDataset) 
484            || anSrcWin[1] + anSrcWin[3] > GDALGetRasterYSize(hDataset) )
485          {
486            fprintf( stderr, 
487                     "Computed -srcwin falls outside raster size of %dx%d.\n",
488                     GDALGetRasterXSize(hDataset), 
489                     GDALGetRasterYSize(hDataset) );
490            exit( 1 );
491          }
492      }
493
494    /* ----------------------------------------------------------------- */
495    /*      Verify source window.                                        */
496    /* ----------------------------------------------------------------- */
497    if( anSrcWin[0] < 0 || anSrcWin[1] < 0 
498        || anSrcWin[2] <= 0 || anSrcWin[3] <= 0
499        || anSrcWin[0] + anSrcWin[2] > GDALGetRasterXSize(hDataset) 
500        || anSrcWin[1] + anSrcWin[3] > GDALGetRasterYSize(hDataset) )
501      {
502        fprintf( stderr, 
503                 "-srcwin %d %d %d %d falls outside raster size of %dx%d\n"
504                 "or is otherwise illegal.\n",
505                 anSrcWin[0],
506                 anSrcWin[1],
507                 anSrcWin[2],
508                 anSrcWin[3],
509                 GDALGetRasterXSize(hDataset), 
510                 GDALGetRasterYSize(hDataset) );
511        exit( 1 );
512      }
513
514    /* ----------------------------------------------------------------- */
515    /*      Find the output driver.                                      */
516    /* ----------------------------------------------------------------- */
517    hDriver = GDALGetDriverByName( pszFormat );
518    if( hDriver == NULL )
519      {
520        int     iDr;
521       
522        char* msg=(char*) CPLMalloc(4096*sizeof(char));
523        sprintf(msg,"Output driver `%s' not recognised.\nThe following format drivers are configured and support output:\n",pszFormat);
524        for( iDr = 0; iDr < GDALGetDriverCount(); iDr++ )
525          {
526            GDALDriverH hDriver = GDALGetDriver(iDr);
527
528            if( GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL ) != NULL
529                || GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATECOPY,
530                                        NULL ) != NULL )
531              {
532                fprintf(stderr,msg);
533                char *tmp=strdup(msg);
534                sprintf( msg,"%s  %s: %s\n",tmp,
535                        GDALGetDriverShortName( hDriver  ),
536                        GDALGetDriverLongName( hDriver ) );
537                free(tmp);
538              }
539          }
540        setMapInMaps(conf,"lenv","message",msg);
541        GDALClose( hDataset );
542        free( panBandList );
543        fflush(stderr);
544        /**
545         * Avoiding GDALDestroyDriverManager() call
546         */
547        CSLDestroy( papszCreateOptions );
548        return 4;
549      }
550
551    /* ----------------------------------------------------------------- */
552    /*   The short form is to CreateCopy().  We use this if the input    */
553    /*   matches the whole dataset.  Eventually we should rewrite        */
554    /*   this entire program to use virtual datasets to construct a      */
555    /*   virtual input source to copy from.                              */
556    /* ----------------------------------------------------------------- */
557    int bSpatialArrangementPreserved = (
558           anSrcWin[0] == 0 && anSrcWin[1] == 0
559        && anSrcWin[2] == GDALGetRasterXSize(hDataset)
560        && anSrcWin[3] == GDALGetRasterYSize(hDataset)
561        && pszOXSize == NULL && pszOYSize == NULL );
562
563    if( eOutputType == GDT_Unknown
564        && !bScale && CSLCount(papszMetadataOptions) == 0 && bDefBands
565        && bSpatialArrangementPreserved
566        && nGCPCount == 0 && !bGotBounds
567        && pszOutputSRS == NULL && !bSetNoData
568        && nRGBExpand == 0)
569      {
570       
571        hOutDS = GDALCreateCopy( hDriver, pszDest, hDataset, 
572                                 FALSE, (char**)0, 
573                                 pfnProgress, NULL );
574
575        if( hOutDS != NULL )
576          GDALClose( hOutDS );
577       
578        GDALClose( hDataset );
579
580        free( panBandList );
581
582        if( !bSubCall )
583          {
584            GDALDumpOpenDatasets( stderr );
585            /**
586             * Avoiding GDALDestroyDriverManager() call
587             */
588          }
589        fprintf(stderr,"==%s %s %s %d==\n",pszFormat,pszSource,pszDest,__LINE__);
590        fflush(stderr);
591
592        CSLDestroy( papszCreateOptions );
593        fprintf(stderr,"==%s %s %s %d==\n",pszFormat,pszSource,pszDest,__LINE__);
594        fflush(stderr);
595        //outputs=(maps*)CPLMalloc(sizeof(maps*));
596        //outputs->name="OutputedPolygon";
597        //outputs->content=createMap("value",(char*)pszDest);
598        //outputs->next=NULL;
599        dumpMaps(outputs);
600        setMapInMaps(outputs,"Result","value",(char*)pszDest);
601        return SERVICE_SUCCEEDED;
602      }
603
604    /* ----------------------------------------------------------------- */
605    /*      Establish some parameters.                                   */
606    /* ----------------------------------------------------------------- */
607    if( pszOXSize == NULL )
608      {
609        nOXSize = anSrcWin[2];
610        nOYSize = anSrcWin[3];
611      }
612    else
613      {
614        nOXSize = (int) ((pszOXSize[strlen(pszOXSize)-1]=='%' 
615                          ? atof(pszOXSize)/100*anSrcWin[2] : atoi(pszOXSize)));
616        nOYSize = (int) ((pszOYSize[strlen(pszOYSize)-1]=='%' 
617                          ? atof(pszOYSize)/100*anSrcWin[3] : atoi(pszOYSize)));
618      }
619
620    /* ================================================================= */
621    /*      Create a virtual dataset.                                    */
622    /* ================================================================= */
623    VRTDataset *poVDS;
624       
625    /* ----------------------------------------------------------------- */
626    /*      Make a virtual clone.                                        */
627    /* ----------------------------------------------------------------- */
628    poVDS = (VRTDataset *) VRTCreate( nOXSize, nOYSize );
629
630    if( nGCPCount == 0 )
631      {
632        if( pszOutputSRS != NULL )
633          {
634            poVDS->SetProjection( pszOutputSRS );
635          }
636        else
637          {
638            pszProjection = GDALGetProjectionRef( hDataset );
639            if( pszProjection != NULL && strlen(pszProjection) > 0 )
640              poVDS->SetProjection( pszProjection );
641          }
642      }
643
644    if( bGotBounds )
645      {
646        adfGeoTransform[0] = adfULLR[0];
647        adfGeoTransform[1] = (adfULLR[2] - adfULLR[0]) / nOXSize;
648        adfGeoTransform[2] = 0.0;
649        adfGeoTransform[3] = adfULLR[1];
650        adfGeoTransform[4] = 0.0;
651        adfGeoTransform[5] = (adfULLR[3] - adfULLR[1]) / nOYSize;
652
653        poVDS->SetGeoTransform( adfGeoTransform );
654      }
655
656    else if( GDALGetGeoTransform( hDataset, adfGeoTransform ) == CE_None
657             && nGCPCount == 0 )
658      {
659        adfGeoTransform[0] += anSrcWin[0] * adfGeoTransform[1]
660          + anSrcWin[1] * adfGeoTransform[2];
661        adfGeoTransform[3] += anSrcWin[0] * adfGeoTransform[4]
662          + anSrcWin[1] * adfGeoTransform[5];
663       
664        adfGeoTransform[1] *= anSrcWin[2] / (double) nOXSize;
665        adfGeoTransform[2] *= anSrcWin[3] / (double) nOYSize;
666        adfGeoTransform[4] *= anSrcWin[2] / (double) nOXSize;
667        adfGeoTransform[5] *= anSrcWin[3] / (double) nOYSize;
668       
669        poVDS->SetGeoTransform( adfGeoTransform );
670      }
671
672    if( nGCPCount != 0 )
673      {
674        const char *pszGCPProjection = pszOutputSRS;
675
676        if( pszGCPProjection == NULL )
677          pszGCPProjection = GDALGetGCPProjection( hDataset );
678        if( pszGCPProjection == NULL )
679          pszGCPProjection = "";
680
681        poVDS->SetGCPs( nGCPCount, pasGCPs, pszGCPProjection );
682
683        GDALDeinitGCPs( nGCPCount, pasGCPs );
684        free( pasGCPs );
685      }
686
687    else if( GDALGetGCPCount( hDataset ) > 0 )
688      {
689        GDAL_GCP *pasGCPs;
690        int       nGCPs = GDALGetGCPCount( hDataset );
691
692        pasGCPs = GDALDuplicateGCPs( nGCPs, GDALGetGCPs( hDataset ) );
693
694        for( i = 0; i < nGCPs; i++ )
695          {
696            pasGCPs[i].dfGCPPixel -= anSrcWin[0];
697            pasGCPs[i].dfGCPLine  -= anSrcWin[1];
698            pasGCPs[i].dfGCPPixel *= (nOXSize / (double) anSrcWin[2] );
699            pasGCPs[i].dfGCPLine  *= (nOYSize / (double) anSrcWin[3] );
700          }
701           
702        poVDS->SetGCPs( nGCPs, pasGCPs,
703                        GDALGetGCPProjection( hDataset ) );
704
705        GDALDeinitGCPs( nGCPs, pasGCPs );
706        free( pasGCPs );
707      }
708
709    /* ----------------------------------------------------------------- */
710    /*      Transfer generally applicable metadata.                      */
711    /* ----------------------------------------------------------------- */
712    poVDS->SetMetadata( ((GDALDataset*)hDataset)->GetMetadata() );
713    AttachMetadata( (GDALDatasetH) poVDS, papszMetadataOptions );
714
715    /* ----------------------------------------------------------------- */
716    /*      Transfer metadata that remains valid if the spatial          */
717    /*      arrangement of the data is unaltered.                        */
718    /* ----------------------------------------------------------------- */
719    if( anSrcWin[0] == 0 && anSrcWin[1] == 0 
720        && anSrcWin[2] == GDALGetRasterXSize(hDataset)
721        && anSrcWin[3] == GDALGetRasterYSize(hDataset) 
722        && pszOXSize == NULL && pszOYSize == NULL )
723      {
724        char **papszMD;
725
726        papszMD = ((GDALDataset*)hDataset)->GetMetadata("RPC");
727        if( papszMD != NULL )
728          poVDS->SetMetadata( papszMD, "RPC" );
729      }
730
731    if (nRGBExpand != 0)
732      nBandCount += nRGBExpand - 1;
733
734    /* ================================================================= */
735    /*      Process all bands.                                           */
736    /* ================================================================= */
737    for( i = 0; i < nBandCount; i++ )
738      {
739        VRTSourcedRasterBand   *poVRTBand;
740        GDALRasterBand  *poSrcBand;
741        GDALDataType    eBandType;
742
743        if (nRGBExpand != 0 && i < nRGBExpand)
744          {
745            poSrcBand = ((GDALDataset *) 
746                         hDataset)->GetRasterBand(panBandList[0]);
747            if (poSrcBand->GetColorTable() == NULL)
748              {
749                fprintf(stderr, "Error : band %d has no color table\n", panBandList[0]);
750                GDALClose( hDataset );
751                free( panBandList );
752                fflush(stderr);
753                /**
754                 * Avoiding GDALDestroyDriverManager() call
755                 */
756                CSLDestroy( papszCreateOptions );
757                exit( 1 );
758              }
759          }
760        else
761          poSrcBand = ((GDALDataset *) 
762                       hDataset)->GetRasterBand(panBandList[i]);
763
764        /* ------------------------------------------------------------ */
765        /*      Select output data type to match source.                */
766        /* ------------------------------------------------------------ */
767        if( eOutputType == GDT_Unknown )
768          eBandType = poSrcBand->GetRasterDataType();
769        else
770          eBandType = eOutputType;
771
772        /* ------------------------------------------------------------ */
773        /*      Create this band.                                       */
774        /* ------------------------------------------------------------ */
775        poVDS->AddBand( eBandType, NULL );
776        poVRTBand = (VRTSourcedRasterBand *) poVDS->GetRasterBand( i+1 );
777           
778        /* ------------------------------------------------------------ */
779        /*      Do we need to collect scaling information?              */
780        /* ------------------------------------------------------------ */
781        double dfScale=1.0, dfOffset=0.0;
782
783        if( bScale && !bHaveScaleSrc )
784          {
785            double      adfCMinMax[2];
786            GDALComputeRasterMinMax( poSrcBand, TRUE, adfCMinMax );
787            dfScaleSrcMin = adfCMinMax[0];
788            dfScaleSrcMax = adfCMinMax[1];
789          }
790
791        if( bScale )
792          {
793            if( dfScaleSrcMax == dfScaleSrcMin )
794              dfScaleSrcMax += 0.1;
795            if( dfScaleDstMax == dfScaleDstMin )
796              dfScaleDstMax += 0.1;
797
798            dfScale = (dfScaleDstMax - dfScaleDstMin) 
799              / (dfScaleSrcMax - dfScaleSrcMin);
800            dfOffset = -1 * dfScaleSrcMin * dfScale + dfScaleDstMin;
801          }
802
803        /* ------------------------------------------------------------ */
804        /*      Create a simple or complex data source depending on the */
805        /*      translation type required.                              */
806        /* ------------------------------------------------------------ */
807        if( bScale || (nRGBExpand != 0 && i < nRGBExpand) )
808          {
809            poVRTBand->AddComplexSource( poSrcBand,
810                                         anSrcWin[0], anSrcWin[1], 
811                                         anSrcWin[2], anSrcWin[3], 
812                                         0, 0, nOXSize, nOYSize,
813                                         dfOffset, dfScale,
814                                         VRT_NODATA_UNSET,
815                                         (nRGBExpand != 0 && i < nRGBExpand) ? i + 1 : 0 );
816          }
817        else
818          poVRTBand->AddSimpleSource( poSrcBand,
819                                      anSrcWin[0], anSrcWin[1], 
820                                      anSrcWin[2], anSrcWin[3], 
821                                      0, 0, nOXSize, nOYSize );
822
823        /* In case of color table translate, we only set the color interpretation */
824        /* other info copied by CopyCommonInfoFrom are not relevant in RGB expansion */
825        if (nRGBExpand != 0 && i < nRGBExpand)
826          {
827            poVRTBand->SetColorInterpretation( (GDALColorInterp) (GCI_RedBand + i) );
828          }
829        else
830          {
831            /* --------------------------------------------------------- */
832            /*      copy over some other information of interest.        */
833            /* --------------------------------------------------------- */
834            poVRTBand->CopyCommonInfoFrom( poSrcBand );
835          }
836
837        /* ------------------------------------------------------------- */
838        /*      Set a forcable nodata value?                             */
839        /* ------------------------------------------------------------- */
840        if( bSetNoData )
841          poVRTBand->SetNoDataValue( dfNoDataReal );
842      }
843
844    /* ----------------------------------------------------------------- */
845    /*      Write to the output file using CopyCreate().                 */
846    /* ----------------------------------------------------------------- */
847    hOutDS = GDALCreateCopy( hDriver, pszDest, (GDALDatasetH) poVDS,
848                             bStrict, papszCreateOptions, 
849                             pfnProgress, NULL );
850
851    if( hOutDS != NULL )
852      {
853        GDALClose( hOutDS );
854      }
855   
856    GDALClose( (GDALDatasetH) poVDS );
857         
858    GDALClose( hDataset );
859   
860    free( panBandList );
861   
862    free( pszOutputSRS );
863
864    if( !bSubCall )
865      {
866        GDALDumpOpenDatasets( stderr );
867        fflush(stderr);
868        /**
869         * Avoiding GDALDestroyDriverManager() call
870         */
871      }
872
873    CSLDestroy( papszCreateOptions );
874   
875    setMapInMaps(outputs,"Result","value",(char*)pszDest);
876
877    return SERVICE_SUCCEEDED;
878  }
879
880
881  /************************************************************************/
882  /*                           AttachMetadata()                           */
883  /************************************************************************/
884
885  static void AttachMetadata( GDALDatasetH hDS, char **papszMetadataOptions )
886
887  {
888    int nCount = CSLCount(papszMetadataOptions);
889    int i;
890
891    for( i = 0; i < nCount; i++ )
892      {
893        char    *pszKey = NULL;
894        const char *pszValue;
895       
896        pszValue = CPLParseNameValue( papszMetadataOptions[i], &pszKey );
897        GDALSetMetadataItem(hDS,pszKey,pszValue,NULL);
898        free( pszKey );
899      }
900
901    CSLDestroy( papszMetadataOptions );
902  }
903
904}
Note: See TracBrowser for help on using the repository browser.

Search

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