source: branches/branch-1.7/zoo-project/zoo-services/gdal/translate/service.c @ 935

Last change on this file since 935 was 448, checked in by djay, 11 years ago

Add support for srcWin in gdal translate service.

  • Property svn:eol-style set to native
  • Property svn:mime-type set to text/x-csrc
File size: 30.4 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 "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
148    tmpMap=NULL;
149    tmpMap=getMapFromMaps(inputs,"SrcWin","value");
150    if(tmpMap!=NULL){
151      char *tmp=tmpMap->value;
152      char *t=strtok(tmp,",");
153      int cnt=0;
154      while(t!=NULL){
155        anSrcWin[cnt] = atoi(t);
156        t=strtok(NULL,",");
157        cnt++;
158      }
159    }   
160
161    tmpMap=NULL;
162    tmpMap=getMapFromMaps(inputs,"ProjWin","value");
163    if(tmpMap!=NULL){
164      char *tmp=tmpMap->value;
165      char *t=strtok(tmp,",");
166      int cnt=0;
167      while(t!=NULL){
168        switch(cnt){
169        case 0:
170          dfULX = atof(t);
171          break;
172        case 1:
173          dfULY = atof(t);
174          break;
175        case 2:
176          dfLRX = atof(t);
177          break;
178        case 3:
179          dfLRY = atof(t);
180          break;
181        }
182        t=strtok(NULL,",");
183        cnt++;
184      }
185    }
186
187    fprintf(stderr,"DEBUG pszDest : %s %d\n",pszDest,__LINE__);
188    tmpMap=NULL;
189    tmpMap=getMapFromMaps(inputs,"GCP","value");
190    if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0){
191      map* length=getMapFromMaps(inputs,"GCP","length");
192      int len=0;
193      if(length){
194        len=atoi(length->value);
195        int i;
196        maps* currentMaps=getMaps(inputs,"GCP");
197        for(i=0;i<len;i++){
198          char* endptr = NULL;
199          /* -gcp pixel line easting northing [elev] */
200         
201          nGCPCount++;
202          pasGCPs = (GDAL_GCP *) 
203            realloc( pasGCPs, sizeof(GDAL_GCP) * nGCPCount );
204          GDALInitGCPs( 1, pasGCPs + nGCPCount - 1 );
205
206          map* currentMap=getMapArray(currentMaps->content,"value",i);
207
208          char* tmpV=strdup(currentMap->value);
209          char *res=strtok(tmpV,",");
210          int j=0;
211          while(res!=NULL){
212            switch(j){
213            case 0:
214              pasGCPs[nGCPCount-1].dfGCPPixel = CPLAtofM(res);
215              break;
216            case 1:
217              pasGCPs[nGCPCount-1].dfGCPLine = CPLAtofM(res);
218              break;
219            case 2:
220              pasGCPs[nGCPCount-1].dfGCPX = CPLAtofM(res);
221              break;
222            case 3:
223              pasGCPs[nGCPCount-1].dfGCPY = CPLAtofM(res);
224              break;
225            case 4:
226              if(res!=NULL && (strtod(res, &endptr) != 0.0 || res[0] == '0'))
227                if (endptr && *endptr == 0)
228                  pasGCPs[nGCPCount-1].dfGCPZ = CPLAtofM(res);
229              break;
230            }
231            res=strtok(NULL,",");
232            j++;
233          }
234        }
235      }else{
236        char* endptr = NULL;
237        /* -gcp pixel line easting northing [elev] */
238       
239        nGCPCount++;
240        pasGCPs = (GDAL_GCP *) 
241          realloc( pasGCPs, sizeof(GDAL_GCP) * nGCPCount );
242        GDALInitGCPs( 1, pasGCPs + nGCPCount - 1 );
243       
244        char* tmpV=strdup(tmpMap->value);
245        char *res=strtok(tmpV,",");
246        int j=0;
247        while(res!=NULL){
248          switch(j){
249          case 0:
250            pasGCPs[nGCPCount-1].dfGCPPixel = CPLAtofM(res);
251            break;
252          case 1:
253            pasGCPs[nGCPCount-1].dfGCPLine = CPLAtofM(res);
254            break;
255          case 2:
256            pasGCPs[nGCPCount-1].dfGCPX = CPLAtofM(res);
257            break;
258          case 3:
259            pasGCPs[nGCPCount-1].dfGCPY = CPLAtofM(res);
260            break;
261          case 4:
262            if(res!=NULL && (CPLStrtod(res, &endptr) != 0.0 || res[0] == '0'))
263              if (endptr && *endptr == 0)
264                pasGCPs[nGCPCount-1].dfGCPZ = CPLAtofM(res);
265            break;
266          }         
267          res=strtok(NULL,",");
268          j++;
269        }
270      }
271    }
272
273    fprintf(stderr,"DEBUG pszDest : %s %d\n",pszDest,__LINE__);
274    tmpMap=NULL;
275    tmpMap=getMapFromMaps(inputs,"SRS","value");
276    if(tmpMap!=NULL){
277      OGRSpatialReference oOutputSRS;
278      if( oOutputSRS.SetFromUserInput( tmpMap->value ) != OGRERR_NONE )
279        {
280          char *msg=(char*)CPLMalloc(100*sizeof(char));
281          sprintf( msg, "Failed to process SRS definition: %s\n", 
282                   tmpMap->value );
283          setMapInMaps(conf,"lenv","message",msg);
284          /**
285           * Avoiding GDALDestroyDriverManager() call
286           */
287          return SERVICE_FAILED;
288        }
289      oOutputSRS.exportToWkt( &pszOutputSRS );
290    }
291    tmpMap=NULL;
292    tmpMap=getMapFromMaps(inputs,"Type","value");
293    if(tmpMap!=NULL){
294      int       iType;
295     
296      for( iType = 1; iType < GDT_TypeCount; iType++ )
297        {
298          if( GDALGetDataTypeName((GDALDataType)iType) != NULL
299              && EQUAL(GDALGetDataTypeName((GDALDataType)iType),
300                       tmpMap->value) )
301            {
302              eOutputType = (GDALDataType) iType;
303            }
304        }
305     
306      if( eOutputType == GDT_Unknown )
307        {
308          printf( "Unknown output pixel type: %s\n", tmpMap->value );
309          /**
310           * Avoiding GDALDestroyDriverManager() call
311           */
312          exit( 2 );
313        }
314    }
315
316    fprintf(stderr,"DEBUG pszDest : %s %d\n",pszDest,__LINE__);
317    if( pszDest == NULL ){
318        fprintf(stderr,"exit line 416");
319        fflush(stderr);
320        /**
321         * Avoiding GDALDestroyDriverManager() call
322         */
323        exit( 10 );
324      }
325
326    fprintf(stderr,"DEBUG pszDest : %s %d\n",pszDest,__LINE__);
327    if ( strcmp(pszSource, pszDest) == 0)
328      {
329        fprintf(stderr, "Source and destination datasets must be different.\n");
330        fflush(stderr);
331        /**
332         * Avoiding GDALDestroyDriverManager() call
333         */
334        exit( 1 );
335      }
336
337    fprintf(stderr,"DEBUG pszDest : %s %d\n",pszDest,__LINE__);
338    /* ----------------------------------------------------------------- */
339    /*      Attempt to open source file.                                 */
340    /* ----------------------------------------------------------------- */
341
342    hDataset = GDALOpenShared( pszSource, GA_ReadOnly );
343    fprintf(stderr,"DEBUG pszDest : %s %d\n",pszDest,__LINE__);
344   
345    if( hDataset == NULL ){
346      char *msg=(char*) CPLMalloc(1024*sizeof(char));
347      sprintf( msg,
348               "GDALOpen failed - %d\n%s\n",
349               CPLGetLastErrorNo(), CPLGetLastErrorMsg() );
350      setMapInMaps(conf,"lenv","message",msg);
351      return SERVICE_FAILED;
352    }
353
354    fprintf(stderr,"DEBUG pszDest : %s %d\n",pszDest,__LINE__);
355    /* ----------------------------------------------------------------- */
356    /*      Handle subdatasets.                                          */
357    /* ----------------------------------------------------------------- */
358    if( !bCopySubDatasets
359        && CSLCount(GDALGetMetadata( hDataset, "SUBDATASETS" )) > 0 
360        && GDALGetRasterCount(hDataset) == 0 )
361      {
362        char *msg=(char*) CPLMalloc(1024*sizeof(char));
363        sprintf( msg,
364                 "Input file contains subdatasets. Please, select one of them for reading.\n" );
365        setMapInMaps(conf,"lenv","message",msg);
366        GDALClose( hDataset );
367        /**
368         * Avoiding GDALDestroyDriverManager() call
369         */
370        exit( 1 );
371      }
372
373    if( CSLCount(GDALGetMetadata( hDataset, "SUBDATASETS" )) > 0 
374        && bCopySubDatasets )
375      {
376        char **papszSubdatasets = GDALGetMetadata(hDataset,"SUBDATASETS");
377        char *pszSubDest = (char *) CPLMalloc(strlen(pszDest)+32);
378        int i;
379        int bOldSubCall = bSubCall;
380       
381        //argv[iDstFileArg] = pszSubDest;
382        bSubCall = TRUE;
383        for( i = 0; papszSubdatasets[i] != NULL; i += 2 )
384          {
385            //argv[iSrcFileArg] = strstr(papszSubdatasets[i],"=")+1;
386            sprintf( pszSubDest, "%s%d", pszDest, i/2 + 1 );
387            /*if( ProxyMain( argc, argv ) != 0 )
388              break;*/
389          }
390       
391        bSubCall = bOldSubCall;
392        free( pszSubDest );
393
394        GDALClose( hDataset );
395
396        if( !bSubCall )
397          {
398            GDALDumpOpenDatasets( stderr );
399            fflush(stderr);
400            /**
401             * Avoiding GDALDestroyDriverManager() call
402             */
403          }
404        return 1;
405      }
406
407    /* ----------------------------------------------------------------- */
408    /*      Collect some information from the source file.               */
409    /* ----------------------------------------------------------------- */
410    nRasterXSize = GDALGetRasterXSize( hDataset );
411    nRasterYSize = GDALGetRasterYSize( hDataset );
412
413    if( !bQuiet )
414      fprintf( stderr, "Input file size is %d, %d\n", nRasterXSize, nRasterYSize );
415
416    if( anSrcWin[2] == 0 && anSrcWin[3] == 0 ){
417        anSrcWin[2] = nRasterXSize;
418        anSrcWin[3] = nRasterYSize;
419      }
420
421    /* ----------------------------------------------------------------- */
422    /*  Build band list to translate                                     */
423    /* ----------------------------------------------------------------- */
424    if( nBandCount == 0 ){
425        nBandCount = GDALGetRasterCount( hDataset );
426        if( nBandCount == 0 ){
427            fprintf( stderr, "Input file has no bands, and so cannot be translated.\n" );
428            fflush(stderr);
429            /**
430             * Avoiding GDALDestroyDriverManager() call
431             */
432            exit(1 );
433          }
434
435        panBandList = (int *) CPLMalloc(sizeof(int)*nBandCount);
436        for( i = 0; i < nBandCount; i++ )
437          panBandList[i] = i+1;
438      }
439    else
440      {
441        for( i = 0; i < nBandCount; i++ )
442          {
443            if( panBandList[i] < 1 || panBandList[i] > GDALGetRasterCount(hDataset) )
444              {
445                fprintf( stderr, 
446                         "Band %d requested, but only bands 1 to %d available.\n",
447                         panBandList[i], GDALGetRasterCount(hDataset) );
448                fflush(stderr);
449                /**
450                 * Avoiding GDALDestroyDriverManager() call
451                 */
452                exit( 2 );
453              }
454          }
455
456        if( nBandCount != GDALGetRasterCount( hDataset ) )
457          bDefBands = FALSE;
458      }
459
460    /* ----------------------------------------------------------------- */
461    /*   Compute the source window from the projected source window      */
462    /*   if the projected coordinates were provided.  Note that the      */
463    /*   projected coordinates are in ulx, uly, lrx, lry format,         */
464    /*   while the anSrcWin is xoff, yoff, xsize, ysize with the         */
465    /*   xoff,yoff being the ulx, uly in pixel/line.                     */
466    /* ----------------------------------------------------------------- */
467    if( dfULX != 0.0 || dfULY != 0.0 
468        || dfLRX != 0.0 || dfLRY != 0.0 )
469      {
470        double  adfGeoTransform[6];
471
472        GDALGetGeoTransform( hDataset, adfGeoTransform );
473
474        if( adfGeoTransform[2] != 0.0 || adfGeoTransform[4] != 0.0 ){
475            fprintf( stderr, 
476                     "The -projwin option was used, but the geotransform is\n"
477                     "rotated.  This configuration is not supported.\n" );
478            GDALClose( hDataset );
479            free( panBandList );
480            fflush(stderr);
481            /**
482             * Avoiding GDALDestroyDriverManager() call
483             */
484            exit( 1 );
485          }
486
487        anSrcWin[0] = (int) 
488          ((dfULX - adfGeoTransform[0]) / adfGeoTransform[1] + 0.001);
489        anSrcWin[1] = (int) 
490          ((dfULY - adfGeoTransform[3]) / adfGeoTransform[5] + 0.001);
491
492        anSrcWin[2] = (int) ((dfLRX - dfULX) / adfGeoTransform[1] + 0.5);
493        anSrcWin[3] = (int) ((dfLRY - dfULY) / adfGeoTransform[5] + 0.5);
494
495        if( !bQuiet )
496          fprintf( stdout, 
497                   "Computed -srcwin %d %d %d %d from projected window.\n",
498                   anSrcWin[0], 
499                   anSrcWin[1], 
500                   anSrcWin[2], 
501                   anSrcWin[3] );
502       
503        if( anSrcWin[0] < 0 || anSrcWin[1] < 0 
504            || anSrcWin[0] + anSrcWin[2] > GDALGetRasterXSize(hDataset) 
505            || anSrcWin[1] + anSrcWin[3] > GDALGetRasterYSize(hDataset) )
506          {
507            fprintf( stderr, 
508                     "Computed -srcwin falls outside raster size of %dx%d.\n",
509                     GDALGetRasterXSize(hDataset), 
510                     GDALGetRasterYSize(hDataset) );
511            exit( 1 );
512          }
513      }
514
515    /* ----------------------------------------------------------------- */
516    /*      Verify source window.                                        */
517    /* ----------------------------------------------------------------- */
518    if( anSrcWin[0] < 0 || anSrcWin[1] < 0 
519        || anSrcWin[2] <= 0 || anSrcWin[3] <= 0
520        || anSrcWin[0] + anSrcWin[2] > GDALGetRasterXSize(hDataset) 
521        || anSrcWin[1] + anSrcWin[3] > GDALGetRasterYSize(hDataset) )
522      {
523        fprintf( stderr, 
524                 "-srcwin %d %d %d %d falls outside raster size of %dx%d\n"
525                 "or is otherwise illegal.\n",
526                 anSrcWin[0],
527                 anSrcWin[1],
528                 anSrcWin[2],
529                 anSrcWin[3],
530                 GDALGetRasterXSize(hDataset), 
531                 GDALGetRasterYSize(hDataset) );
532        exit( 1 );
533      }
534
535    /* ----------------------------------------------------------------- */
536    /*      Find the output driver.                                      */
537    /* ----------------------------------------------------------------- */
538    hDriver = GDALGetDriverByName( pszFormat );
539    if( hDriver == NULL )
540      {
541        int     iDr;
542       
543        char* msg=(char*) CPLMalloc(4096*sizeof(char));
544        sprintf(msg,"Output driver `%s' not recognised.\nThe following format drivers are configured and support output:\n",pszFormat);
545        for( iDr = 0; iDr < GDALGetDriverCount(); iDr++ )
546          {
547            GDALDriverH hDriver = GDALGetDriver(iDr);
548
549            if( GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL ) != NULL
550                || GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATECOPY,
551                                        NULL ) != NULL )
552              {
553                fprintf(stderr,msg);
554                char *tmp=strdup(msg);
555                sprintf( msg,"%s  %s: %s\n",tmp,
556                        GDALGetDriverShortName( hDriver  ),
557                        GDALGetDriverLongName( hDriver ) );
558                free(tmp);
559              }
560          }
561        setMapInMaps(conf,"lenv","message",msg);
562        GDALClose( hDataset );
563        free( panBandList );
564        fflush(stderr);
565        /**
566         * Avoiding GDALDestroyDriverManager() call
567         */
568        CSLDestroy( papszCreateOptions );
569        return 4;
570      }
571
572    /* ----------------------------------------------------------------- */
573    /*   The short form is to CreateCopy().  We use this if the input    */
574    /*   matches the whole dataset.  Eventually we should rewrite        */
575    /*   this entire program to use virtual datasets to construct a      */
576    /*   virtual input source to copy from.                              */
577    /* ----------------------------------------------------------------- */
578    int bSpatialArrangementPreserved = (
579           anSrcWin[0] == 0 && anSrcWin[1] == 0
580        && anSrcWin[2] == GDALGetRasterXSize(hDataset)
581        && anSrcWin[3] == GDALGetRasterYSize(hDataset)
582        && pszOXSize == NULL && pszOYSize == NULL );
583
584    if( eOutputType == GDT_Unknown
585        && !bScale && CSLCount(papszMetadataOptions) == 0 && bDefBands
586        && bSpatialArrangementPreserved
587        && nGCPCount == 0 && !bGotBounds
588        && pszOutputSRS == NULL && !bSetNoData
589        && nRGBExpand == 0)
590      {
591       
592        hOutDS = GDALCreateCopy( hDriver, pszDest, hDataset, 
593                                 FALSE, (char**)0, 
594                                 pfnProgress, NULL );
595
596        if( hOutDS != NULL )
597          GDALClose( hOutDS );
598       
599        GDALClose( hDataset );
600
601        free( panBandList );
602
603        if( !bSubCall )
604          {
605            GDALDumpOpenDatasets( stderr );
606            /**
607             * Avoiding GDALDestroyDriverManager() call
608             */
609          }
610        fprintf(stderr,"==%s %s %s %d==\n",pszFormat,pszSource,pszDest,__LINE__);
611        fflush(stderr);
612
613        CSLDestroy( papszCreateOptions );
614        fprintf(stderr,"==%s %s %s %d==\n",pszFormat,pszSource,pszDest,__LINE__);
615        fflush(stderr);
616        //outputs=(maps*)CPLMalloc(sizeof(maps*));
617        //outputs->name="OutputedPolygon";
618        //outputs->content=createMap("value",(char*)pszDest);
619        //outputs->next=NULL;
620        dumpMaps(outputs);
621        setMapInMaps(outputs,"Result","value",(char*)pszDest);
622        return SERVICE_SUCCEEDED;
623      }
624
625    /* ----------------------------------------------------------------- */
626    /*      Establish some parameters.                                   */
627    /* ----------------------------------------------------------------- */
628    if( pszOXSize == NULL )
629      {
630        nOXSize = anSrcWin[2];
631        nOYSize = anSrcWin[3];
632      }
633    else
634      {
635        nOXSize = (int) ((pszOXSize[strlen(pszOXSize)-1]=='%' 
636                          ? atof(pszOXSize)/100*anSrcWin[2] : atoi(pszOXSize)));
637        nOYSize = (int) ((pszOYSize[strlen(pszOYSize)-1]=='%' 
638                          ? atof(pszOYSize)/100*anSrcWin[3] : atoi(pszOYSize)));
639      }
640
641    /* ================================================================= */
642    /*      Create a virtual dataset.                                    */
643    /* ================================================================= */
644    VRTDataset *poVDS;
645       
646    /* ----------------------------------------------------------------- */
647    /*      Make a virtual clone.                                        */
648    /* ----------------------------------------------------------------- */
649    poVDS = (VRTDataset *) VRTCreate( nOXSize, nOYSize );
650
651    if( nGCPCount == 0 )
652      {
653        if( pszOutputSRS != NULL )
654          {
655            poVDS->SetProjection( pszOutputSRS );
656          }
657        else
658          {
659            pszProjection = GDALGetProjectionRef( hDataset );
660            if( pszProjection != NULL && strlen(pszProjection) > 0 )
661              poVDS->SetProjection( pszProjection );
662          }
663      }
664
665    if( bGotBounds )
666      {
667        adfGeoTransform[0] = adfULLR[0];
668        adfGeoTransform[1] = (adfULLR[2] - adfULLR[0]) / nOXSize;
669        adfGeoTransform[2] = 0.0;
670        adfGeoTransform[3] = adfULLR[1];
671        adfGeoTransform[4] = 0.0;
672        adfGeoTransform[5] = (adfULLR[3] - adfULLR[1]) / nOYSize;
673
674        poVDS->SetGeoTransform( adfGeoTransform );
675      }
676
677    else if( GDALGetGeoTransform( hDataset, adfGeoTransform ) == CE_None
678             && nGCPCount == 0 )
679      {
680        adfGeoTransform[0] += anSrcWin[0] * adfGeoTransform[1]
681          + anSrcWin[1] * adfGeoTransform[2];
682        adfGeoTransform[3] += anSrcWin[0] * adfGeoTransform[4]
683          + anSrcWin[1] * adfGeoTransform[5];
684       
685        adfGeoTransform[1] *= anSrcWin[2] / (double) nOXSize;
686        adfGeoTransform[2] *= anSrcWin[3] / (double) nOYSize;
687        adfGeoTransform[4] *= anSrcWin[2] / (double) nOXSize;
688        adfGeoTransform[5] *= anSrcWin[3] / (double) nOYSize;
689       
690        poVDS->SetGeoTransform( adfGeoTransform );
691      }
692
693    if( nGCPCount != 0 )
694      {
695        const char *pszGCPProjection = pszOutputSRS;
696
697        if( pszGCPProjection == NULL )
698          pszGCPProjection = GDALGetGCPProjection( hDataset );
699        if( pszGCPProjection == NULL )
700          pszGCPProjection = "";
701
702        poVDS->SetGCPs( nGCPCount, pasGCPs, pszGCPProjection );
703
704        GDALDeinitGCPs( nGCPCount, pasGCPs );
705        free( pasGCPs );
706      }
707
708    else if( GDALGetGCPCount( hDataset ) > 0 )
709      {
710        GDAL_GCP *pasGCPs;
711        int       nGCPs = GDALGetGCPCount( hDataset );
712
713        pasGCPs = GDALDuplicateGCPs( nGCPs, GDALGetGCPs( hDataset ) );
714
715        for( i = 0; i < nGCPs; i++ )
716          {
717            pasGCPs[i].dfGCPPixel -= anSrcWin[0];
718            pasGCPs[i].dfGCPLine  -= anSrcWin[1];
719            pasGCPs[i].dfGCPPixel *= (nOXSize / (double) anSrcWin[2] );
720            pasGCPs[i].dfGCPLine  *= (nOYSize / (double) anSrcWin[3] );
721          }
722           
723        poVDS->SetGCPs( nGCPs, pasGCPs,
724                        GDALGetGCPProjection( hDataset ) );
725
726        GDALDeinitGCPs( nGCPs, pasGCPs );
727        free( pasGCPs );
728      }
729
730    /* ----------------------------------------------------------------- */
731    /*      Transfer generally applicable metadata.                      */
732    /* ----------------------------------------------------------------- */
733    poVDS->SetMetadata( ((GDALDataset*)hDataset)->GetMetadata() );
734    AttachMetadata( (GDALDatasetH) poVDS, papszMetadataOptions );
735
736    /* ----------------------------------------------------------------- */
737    /*      Transfer metadata that remains valid if the spatial          */
738    /*      arrangement of the data is unaltered.                        */
739    /* ----------------------------------------------------------------- */
740    if( anSrcWin[0] == 0 && anSrcWin[1] == 0 
741        && anSrcWin[2] == GDALGetRasterXSize(hDataset)
742        && anSrcWin[3] == GDALGetRasterYSize(hDataset) 
743        && pszOXSize == NULL && pszOYSize == NULL )
744      {
745        char **papszMD;
746
747        papszMD = ((GDALDataset*)hDataset)->GetMetadata("RPC");
748        if( papszMD != NULL )
749          poVDS->SetMetadata( papszMD, "RPC" );
750      }
751
752    if (nRGBExpand != 0)
753      nBandCount += nRGBExpand - 1;
754
755    /* ================================================================= */
756    /*      Process all bands.                                           */
757    /* ================================================================= */
758    for( i = 0; i < nBandCount; i++ )
759      {
760        VRTSourcedRasterBand   *poVRTBand;
761        GDALRasterBand  *poSrcBand;
762        GDALDataType    eBandType;
763
764        if (nRGBExpand != 0 && i < nRGBExpand)
765          {
766            poSrcBand = ((GDALDataset *) 
767                         hDataset)->GetRasterBand(panBandList[0]);
768            if (poSrcBand->GetColorTable() == NULL)
769              {
770                fprintf(stderr, "Error : band %d has no color table\n", panBandList[0]);
771                GDALClose( hDataset );
772                free( panBandList );
773                fflush(stderr);
774                /**
775                 * Avoiding GDALDestroyDriverManager() call
776                 */
777                CSLDestroy( papszCreateOptions );
778                exit( 1 );
779              }
780          }
781        else
782          poSrcBand = ((GDALDataset *) 
783                       hDataset)->GetRasterBand(panBandList[i]);
784
785        /* ------------------------------------------------------------ */
786        /*      Select output data type to match source.                */
787        /* ------------------------------------------------------------ */
788        if( eOutputType == GDT_Unknown )
789          eBandType = poSrcBand->GetRasterDataType();
790        else
791          eBandType = eOutputType;
792
793        /* ------------------------------------------------------------ */
794        /*      Create this band.                                       */
795        /* ------------------------------------------------------------ */
796        poVDS->AddBand( eBandType, NULL );
797        poVRTBand = (VRTSourcedRasterBand *) poVDS->GetRasterBand( i+1 );
798           
799        /* ------------------------------------------------------------ */
800        /*      Do we need to collect scaling information?              */
801        /* ------------------------------------------------------------ */
802        double dfScale=1.0, dfOffset=0.0;
803
804        if( bScale && !bHaveScaleSrc )
805          {
806            double      adfCMinMax[2];
807            GDALComputeRasterMinMax( poSrcBand, TRUE, adfCMinMax );
808            dfScaleSrcMin = adfCMinMax[0];
809            dfScaleSrcMax = adfCMinMax[1];
810          }
811
812        if( bScale )
813          {
814            if( dfScaleSrcMax == dfScaleSrcMin )
815              dfScaleSrcMax += 0.1;
816            if( dfScaleDstMax == dfScaleDstMin )
817              dfScaleDstMax += 0.1;
818
819            dfScale = (dfScaleDstMax - dfScaleDstMin) 
820              / (dfScaleSrcMax - dfScaleSrcMin);
821            dfOffset = -1 * dfScaleSrcMin * dfScale + dfScaleDstMin;
822          }
823
824        /* ------------------------------------------------------------ */
825        /*      Create a simple or complex data source depending on the */
826        /*      translation type required.                              */
827        /* ------------------------------------------------------------ */
828        if( bScale || (nRGBExpand != 0 && i < nRGBExpand) )
829          {
830            poVRTBand->AddComplexSource( poSrcBand,
831                                         anSrcWin[0], anSrcWin[1], 
832                                         anSrcWin[2], anSrcWin[3], 
833                                         0, 0, nOXSize, nOYSize,
834                                         dfOffset, dfScale,
835                                         VRT_NODATA_UNSET,
836                                         (nRGBExpand != 0 && i < nRGBExpand) ? i + 1 : 0 );
837          }
838        else
839          poVRTBand->AddSimpleSource( poSrcBand,
840                                      anSrcWin[0], anSrcWin[1], 
841                                      anSrcWin[2], anSrcWin[3], 
842                                      0, 0, nOXSize, nOYSize );
843
844        /* In case of color table translate, we only set the color interpretation */
845        /* other info copied by CopyCommonInfoFrom are not relevant in RGB expansion */
846        if (nRGBExpand != 0 && i < nRGBExpand)
847          {
848            poVRTBand->SetColorInterpretation( (GDALColorInterp) (GCI_RedBand + i) );
849          }
850        else
851          {
852            /* --------------------------------------------------------- */
853            /*      copy over some other information of interest.        */
854            /* --------------------------------------------------------- */
855            poVRTBand->CopyCommonInfoFrom( poSrcBand );
856          }
857
858        /* ------------------------------------------------------------- */
859        /*      Set a forcable nodata value?                             */
860        /* ------------------------------------------------------------- */
861        if( bSetNoData )
862          poVRTBand->SetNoDataValue( dfNoDataReal );
863      }
864
865    /* ----------------------------------------------------------------- */
866    /*      Write to the output file using CopyCreate().                 */
867    /* ----------------------------------------------------------------- */
868    hOutDS = GDALCreateCopy( hDriver, pszDest, (GDALDatasetH) poVDS,
869                             bStrict, papszCreateOptions, 
870                             pfnProgress, NULL );
871
872    if( hOutDS != NULL )
873      {
874        GDALClose( hOutDS );
875      }
876   
877    GDALClose( (GDALDatasetH) poVDS );
878         
879    GDALClose( hDataset );
880   
881    free( panBandList );
882   
883    free( pszOutputSRS );
884
885    if( !bSubCall )
886      {
887        GDALDumpOpenDatasets( stderr );
888        fflush(stderr);
889        /**
890         * Avoiding GDALDestroyDriverManager() call
891         */
892      }
893
894    CSLDestroy( papszCreateOptions );
895   
896    setMapInMaps(outputs,"Result","value",(char*)pszDest);
897
898    return SERVICE_SUCCEEDED;
899  }
900
901
902  /************************************************************************/
903  /*                           AttachMetadata()                           */
904  /************************************************************************/
905
906  static void AttachMetadata( GDALDatasetH hDS, char **papszMetadataOptions )
907
908  {
909    int nCount = CSLCount(papszMetadataOptions);
910    int i;
911
912    for( i = 0; i < nCount; i++ )
913      {
914        char    *pszKey = NULL;
915        const char *pszValue;
916       
917        pszValue = CPLParseNameValue( papszMetadataOptions[i], &pszKey );
918        GDALSetMetadataItem(hDS,pszKey,pszValue,NULL);
919        free( pszKey );
920      }
921
922    CSLDestroy( papszMetadataOptions );
923  }
924
925}
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