1 : /******************************************************************************
2 : * $Id: vrtwarped.cpp 18949 2010-02-27 22:54:47Z rouault $
3 : *
4 : * Project: Virtual GDAL Datasets
5 : * Purpose: Implementation of VRTWarpedRasterBand *and VRTWarpedDataset.
6 : * Author: Frank Warmerdam <warmerdam@pobox.com>
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2004, Frank Warmerdam <warmerdam@pobox.com>
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 "vrtdataset.h"
31 : #include "cpl_minixml.h"
32 : #include "cpl_string.h"
33 : #include "gdalwarper.h"
34 : #include <cassert>
35 :
36 : CPL_CVSID("$Id: vrtwarped.cpp 18949 2010-02-27 22:54:47Z rouault $");
37 :
38 : /************************************************************************/
39 : /* GDALAutoCreateWarpedVRT() */
40 : /************************************************************************/
41 :
42 : /**
43 : * Create virtual warped dataset automatically.
44 : *
45 : * This function will create a warped virtual file representing the
46 : * input image warped into the target coordinate system. A GenImgProj
47 : * transformation is created to accomplish any required GCP/Geotransform
48 : * warp and reprojection to the target coordinate system. The output virtual
49 : * dataset will be "northup" in the target coordinate system. The
50 : * GDALSuggestedWarpOutput() function is used to determine the bounds and
51 : * resolution of the output virtual file which should be large enough to
52 : * include all the input image
53 : *
54 : * Note that the constructed GDALDatasetH will acquire one or more references
55 : * to the passed in hSrcDS. Reference counting semantics on the source
56 : * dataset should be honoured. That is, don't just GDALClose() it unless it
57 : * was opened with GDALOpenShared().
58 : *
59 : * The returned dataset will have no associated filename for itself. If you
60 : * want to write the virtual dataset description to a file, use the
61 : * GDALSetDescription() function (or SetDescription() method) on the dataset
62 : * to assign a filename before it is closed.
63 : *
64 : * @param hSrcDS The source dataset.
65 : *
66 : * @param pszSrcWKT The coordinate system of the source image. If NULL, it
67 : * will be read from the source image.
68 : *
69 : * @param pszDstWKT The coordinate system to convert to. If NULL no change
70 : * of coordinate system will take place.
71 : *
72 : * @param eResampleAlg One of GRA_NearestNeighbour, GRA_Bilinear, GRA_Cubic or
73 : * GRA_CubicSpline. Controls the sampling method used.
74 : *
75 : * @param dfMaxError Maximum error measured in input pixels that is allowed in
76 : * approximating the transformation (0.0 for exact calculations).
77 : *
78 : * @param psOptions Additional warp options, normally NULL.
79 : *
80 : * @return NULL on failure, or a new virtual dataset handle on success.
81 : */
82 :
83 : GDALDatasetH CPL_STDCALL
84 : GDALAutoCreateWarpedVRT( GDALDatasetH hSrcDS,
85 : const char *pszSrcWKT,
86 : const char *pszDstWKT,
87 : GDALResampleAlg eResampleAlg,
88 : double dfMaxError,
89 1 : const GDALWarpOptions *psOptionsIn )
90 :
91 : {
92 : GDALWarpOptions *psWO;
93 : int i;
94 :
95 1 : VALIDATE_POINTER1( hSrcDS, "GDALAutoCreateWarpedVRT", NULL );
96 :
97 : /* -------------------------------------------------------------------- */
98 : /* Populate the warp options. */
99 : /* -------------------------------------------------------------------- */
100 1 : if( psOptionsIn != NULL )
101 0 : psWO = GDALCloneWarpOptions( psOptionsIn );
102 : else
103 1 : psWO = GDALCreateWarpOptions();
104 :
105 1 : psWO->eResampleAlg = eResampleAlg;
106 :
107 1 : psWO->hSrcDS = hSrcDS;
108 :
109 1 : psWO->nBandCount = GDALGetRasterCount( hSrcDS );
110 1 : psWO->panSrcBands = (int *) CPLMalloc(sizeof(int) * psWO->nBandCount);
111 1 : psWO->panDstBands = (int *) CPLMalloc(sizeof(int) * psWO->nBandCount);
112 :
113 4 : for( i = 0; i < psWO->nBandCount; i++ )
114 : {
115 3 : psWO->panSrcBands[i] = i+1;
116 3 : psWO->panDstBands[i] = i+1;
117 : }
118 :
119 : /* TODO: should fill in no data where available */
120 :
121 : /* -------------------------------------------------------------------- */
122 : /* Create the transformer. */
123 : /* -------------------------------------------------------------------- */
124 1 : psWO->pfnTransformer = GDALGenImgProjTransform;
125 : psWO->pTransformerArg =
126 : GDALCreateGenImgProjTransformer( psWO->hSrcDS, pszSrcWKT,
127 : NULL, pszDstWKT,
128 1 : TRUE, 1.0, 0 );
129 :
130 1 : if( psWO->pTransformerArg == NULL )
131 : {
132 0 : GDALDestroyWarpOptions( psWO );
133 0 : return NULL;
134 : }
135 :
136 : /* -------------------------------------------------------------------- */
137 : /* Figure out the desired output bounds and resolution. */
138 : /* -------------------------------------------------------------------- */
139 : double adfDstGeoTransform[6];
140 : int nDstPixels, nDstLines;
141 : CPLErr eErr;
142 :
143 : eErr =
144 : GDALSuggestedWarpOutput( hSrcDS, psWO->pfnTransformer,
145 : psWO->pTransformerArg,
146 1 : adfDstGeoTransform, &nDstPixels, &nDstLines );
147 :
148 : /* -------------------------------------------------------------------- */
149 : /* Update the transformer to include an output geotransform */
150 : /* back to pixel/line coordinates. */
151 : /* */
152 : /* -------------------------------------------------------------------- */
153 : GDALSetGenImgProjTransformerDstGeoTransform(
154 1 : psWO->pTransformerArg, adfDstGeoTransform );
155 :
156 : /* -------------------------------------------------------------------- */
157 : /* Do we want to apply an approximating transformation? */
158 : /* -------------------------------------------------------------------- */
159 1 : if( dfMaxError > 0.0 )
160 : {
161 : psWO->pTransformerArg =
162 : GDALCreateApproxTransformer( psWO->pfnTransformer,
163 : psWO->pTransformerArg,
164 0 : dfMaxError );
165 0 : psWO->pfnTransformer = GDALApproxTransform;
166 : }
167 :
168 : /* -------------------------------------------------------------------- */
169 : /* Create the VRT file. */
170 : /* -------------------------------------------------------------------- */
171 : GDALDatasetH hDstDS;
172 :
173 : hDstDS = GDALCreateWarpedVRT( hSrcDS, nDstPixels, nDstLines,
174 1 : adfDstGeoTransform, psWO );
175 :
176 1 : GDALDestroyWarpOptions( psWO );
177 :
178 1 : if( pszDstWKT != NULL )
179 0 : GDALSetProjection( hDstDS, pszDstWKT );
180 1 : else if( pszSrcWKT != NULL )
181 0 : GDALSetProjection( hDstDS, pszDstWKT );
182 1 : else if( GDALGetGCPCount( hSrcDS ) > 0 )
183 1 : GDALSetProjection( hDstDS, GDALGetGCPProjection( hSrcDS ) );
184 : else
185 0 : GDALSetProjection( hDstDS, GDALGetProjectionRef( hSrcDS ) );
186 :
187 1 : return hDstDS;
188 : }
189 :
190 : /************************************************************************/
191 : /* GDALCreateWarpedVRT() */
192 : /************************************************************************/
193 :
194 : /**
195 : * Create virtual warped dataset.
196 : *
197 : * This function will create a warped virtual file representing the
198 : * input image warped based on a provided transformation. Output bounds
199 : * and resolution are provided explicitly.
200 : *
201 : * Note that the constructed GDALDatasetH will acquire one or more references
202 : * to the passed in hSrcDS. Reference counting semantics on the source
203 : * dataset should be honoured. That is, don't just GDALClose() it unless it
204 : * was opened with GDALOpenShared().
205 : *
206 : * @param hSrcDS The source dataset.
207 : *
208 : * @param nPixels Width of the virtual warped dataset to create
209 : *
210 : * @param nLines Height of the virtual warped dataset to create
211 : *
212 : * @param padfGeoTransform Geotransform matrix of the virtual warped dataset to create
213 : *
214 : * @param psOptions Warp options. Must be different from NULL.
215 : *
216 : * @return NULL on failure, or a new virtual dataset handle on success.
217 : */
218 :
219 : GDALDatasetH CPL_STDCALL
220 : GDALCreateWarpedVRT( GDALDatasetH hSrcDS,
221 : int nPixels, int nLines, double *padfGeoTransform,
222 1 : GDALWarpOptions *psOptions )
223 :
224 : {
225 1 : VALIDATE_POINTER1( hSrcDS, "GDALCreateWarpedVRT", NULL );
226 :
227 : /* -------------------------------------------------------------------- */
228 : /* Create the VRTDataset and populate it with bands. */
229 : /* -------------------------------------------------------------------- */
230 1 : VRTWarpedDataset *poDS = new VRTWarpedDataset( nPixels, nLines );
231 : int i;
232 :
233 1 : psOptions->hDstDS = (GDALDatasetH) poDS;
234 :
235 1 : poDS->SetGeoTransform( padfGeoTransform );
236 :
237 4 : for( i = 0; i < psOptions->nBandCount; i++ )
238 : {
239 : VRTWarpedRasterBand *poBand;
240 : GDALRasterBand *poSrcBand = (GDALRasterBand *)
241 3 : GDALGetRasterBand( hSrcDS, i+1 );
242 :
243 3 : poDS->AddBand( poSrcBand->GetRasterDataType(), NULL );
244 :
245 3 : poBand = (VRTWarpedRasterBand *) poDS->GetRasterBand( i+1 );
246 3 : poBand->CopyCommonInfoFrom( poSrcBand );
247 : }
248 :
249 : /* -------------------------------------------------------------------- */
250 : /* Initialize the warp on the VRTWarpedDataset. */
251 : /* -------------------------------------------------------------------- */
252 1 : poDS->Initialize( psOptions );
253 :
254 1 : return (GDALDatasetH) poDS;
255 : }
256 :
257 : /************************************************************************/
258 : /* ==================================================================== */
259 : /* VRTWarpedDataset */
260 : /* ==================================================================== */
261 : /************************************************************************/
262 :
263 : /************************************************************************/
264 : /* VRTWarpedDataset() */
265 : /************************************************************************/
266 :
267 43 : VRTWarpedDataset::VRTWarpedDataset( int nXSize, int nYSize )
268 43 : : VRTDataset( nXSize, nYSize )
269 :
270 : {
271 43 : poWarper = NULL;
272 43 : nBlockXSize = 512;
273 43 : nBlockYSize = 128;
274 43 : eAccess = GA_Update;
275 :
276 43 : nOverviewCount = 0;
277 43 : papoOverviews = NULL;
278 43 : }
279 :
280 : /************************************************************************/
281 : /* ~VRTWarpedDataset() */
282 : /************************************************************************/
283 :
284 43 : VRTWarpedDataset::~VRTWarpedDataset()
285 :
286 : {
287 43 : FlushCache();
288 :
289 : /* -------------------------------------------------------------------- */
290 : /* Cleanup overviews. */
291 : /* -------------------------------------------------------------------- */
292 : int iOverview;
293 :
294 45 : for( iOverview = 0; iOverview < nOverviewCount; iOverview++ )
295 : {
296 2 : GDALDatasetH hDS = (GDALDatasetH) papoOverviews[iOverview];
297 :
298 2 : if( GDALDereferenceDataset( hDS ) < 1 )
299 : {
300 2 : GDALReferenceDataset( hDS );
301 2 : GDALClose( hDS );
302 : }
303 : }
304 :
305 43 : CPLFree( papoOverviews );
306 :
307 : /* -------------------------------------------------------------------- */
308 : /* Cleanup warper if one is in effect. */
309 : /* -------------------------------------------------------------------- */
310 43 : if( poWarper != NULL )
311 : {
312 43 : const GDALWarpOptions *psWO = poWarper->GetOptions();
313 :
314 : /* -------------------------------------------------------------------- */
315 : /* We take care to only call GDALClose() on psWO->hSrcDS if the */
316 : /* reference count drops to zero. This is makes it so that we */
317 : /* can operate reference counting semantics more-or-less */
318 : /* properly even if the dataset isn't open in shared mode, */
319 : /* though we require that the caller also honour the reference */
320 : /* counting semantics even though it isn't a shared dataset. */
321 : /* -------------------------------------------------------------------- */
322 43 : if( psWO->hSrcDS != NULL )
323 : {
324 43 : if( GDALDereferenceDataset( psWO->hSrcDS ) < 1 )
325 : {
326 39 : GDALReferenceDataset( psWO->hSrcDS );
327 39 : GDALClose( psWO->hSrcDS );
328 : }
329 : }
330 :
331 : /* -------------------------------------------------------------------- */
332 : /* We are responsible for cleaning up the transformer outselves. */
333 : /* -------------------------------------------------------------------- */
334 43 : if( psWO->pTransformerArg != NULL )
335 43 : GDALDestroyTransformer( psWO->pTransformerArg );
336 :
337 43 : delete poWarper;
338 : }
339 43 : }
340 :
341 : /************************************************************************/
342 : /* Initialize() */
343 : /* */
344 : /* Initialize a dataset from passed in warp options. */
345 : /************************************************************************/
346 :
347 5 : CPLErr VRTWarpedDataset::Initialize( void *psWO )
348 :
349 : {
350 5 : if( poWarper != NULL )
351 0 : delete poWarper;
352 :
353 5 : poWarper = new GDALWarpOperation();
354 :
355 : // The act of initializing this warped dataset with this warp options
356 : // will result in our assuming ownership of a reference to the
357 : // hSrcDS.
358 :
359 5 : if( ((GDALWarpOptions *) psWO)->hSrcDS != NULL )
360 5 : GDALReferenceDataset( ((GDALWarpOptions *) psWO)->hSrcDS );
361 :
362 5 : return poWarper->Initialize( (GDALWarpOptions *) psWO );
363 : }
364 :
365 : /************************************************************************/
366 : /* GetFileList() */
367 : /************************************************************************/
368 :
369 0 : char** VRTWarpedDataset::GetFileList()
370 : {
371 0 : char** papszFileList = GDALDataset::GetFileList();
372 :
373 0 : if( poWarper != NULL )
374 : {
375 0 : const GDALWarpOptions *psWO = poWarper->GetOptions();
376 : const char* pszFilename;
377 :
378 0 : if( psWO->hSrcDS != NULL &&
379 : (pszFilename =
380 : ((GDALDataset*)psWO->hSrcDS)->GetDescription()) != NULL )
381 : {
382 : VSIStatBufL sStat;
383 0 : if( VSIStatL( pszFilename, &sStat ) == 0 )
384 : {
385 0 : papszFileList = CSLAddString(papszFileList, pszFilename);
386 : }
387 : }
388 : }
389 :
390 0 : return papszFileList;
391 : }
392 :
393 : /************************************************************************/
394 : /* VRTWarpedOverviewTransform() */
395 : /************************************************************************/
396 :
397 : typedef struct {
398 : GDALTransformerInfo sTI;
399 :
400 : GDALTransformerFunc pfnBaseTransformer;
401 : void *pBaseTransformerArg;
402 :
403 : double dfXOverviewFactor;
404 : double dfYOverviewFactor;
405 : } VWOTInfo;
406 :
407 : static
408 : int VRTWarpedOverviewTransform( void *pTransformArg, int bDstToSrc,
409 : int nPointCount,
410 : double *padfX, double *padfY, double *padfZ,
411 258 : int *panSuccess )
412 :
413 : {
414 258 : VWOTInfo *psInfo = (VWOTInfo *) pTransformArg;
415 : int i, bSuccess;
416 :
417 258 : if( bDstToSrc )
418 : {
419 131498 : for( i = 0; i < nPointCount; i++ )
420 : {
421 131240 : padfX[i] *= psInfo->dfXOverviewFactor;
422 131240 : padfY[i] *= psInfo->dfYOverviewFactor;
423 : }
424 : }
425 :
426 : bSuccess = psInfo->pfnBaseTransformer( psInfo->pBaseTransformerArg,
427 : bDstToSrc,
428 : nPointCount, padfX, padfY, padfZ,
429 258 : panSuccess );
430 :
431 258 : if( !bDstToSrc )
432 : {
433 0 : for( i = 0; i < nPointCount; i++ )
434 : {
435 0 : padfX[i] /= psInfo->dfXOverviewFactor;
436 0 : padfY[i] /= psInfo->dfYOverviewFactor;
437 : }
438 : }
439 :
440 258 : return bSuccess;
441 : }
442 :
443 2 : static void VRTWarpedOverviewCleanup(void* pTransformArg)
444 : {
445 2 : VWOTInfo *psInfo = (VWOTInfo *) pTransformArg;
446 :
447 2 : CPLFree( psInfo );
448 2 : }
449 :
450 : /************************************************************************/
451 : /* BuildOverviews() */
452 : /* */
453 : /* For overviews, we actually just build a whole new dataset */
454 : /* with an extra layer of transformation on the warper used to */
455 : /* accomplish downsampling by the desired factor. */
456 : /************************************************************************/
457 :
458 : CPLErr
459 : VRTWarpedDataset::IBuildOverviews( const char *pszResampling,
460 : int nOverviews, int *panOverviewList,
461 : int nListBands, int *panBandList,
462 : GDALProgressFunc pfnProgress,
463 2 : void * pProgressData )
464 :
465 : {
466 : /* -------------------------------------------------------------------- */
467 : /* Initial progress result. */
468 : /* -------------------------------------------------------------------- */
469 2 : if( !pfnProgress( 0.0, NULL, pProgressData ) )
470 : {
471 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
472 0 : return CE_Failure;
473 : }
474 :
475 : /* -------------------------------------------------------------------- */
476 : /* Establish which of the overview levels we already have, and */
477 : /* which are new. */
478 : /* -------------------------------------------------------------------- */
479 2 : int i, nNewOverviews, *panNewOverviewList = NULL;
480 :
481 2 : nNewOverviews = 0;
482 2 : panNewOverviewList = (int *) CPLCalloc(sizeof(int),nOverviews);
483 4 : for( i = 0; i < nOverviews; i++ )
484 : {
485 : int j;
486 :
487 2 : for( j = 0; j < nOverviewCount; j++ )
488 : {
489 : int nOvFactor;
490 0 : VRTWarpedDataset *poOverview = papoOverviews[j];
491 :
492 : nOvFactor = (int)
493 0 : (0.5+GetRasterXSize() / (double) poOverview->GetRasterXSize());
494 :
495 0 : if( nOvFactor == panOverviewList[i]
496 : || nOvFactor == GDALOvLevelAdjust( panOverviewList[i],
497 : GetRasterXSize() ) )
498 0 : panOverviewList[i] *= -1;
499 : }
500 :
501 2 : if( panOverviewList[i] > 0 )
502 2 : panNewOverviewList[nNewOverviews++] = panOverviewList[i];
503 : }
504 :
505 : /* -------------------------------------------------------------------- */
506 : /* Create each missing overview (we don't need to do anything */
507 : /* to update existing overviews). */
508 : /* -------------------------------------------------------------------- */
509 4 : for( i = 0; i < nNewOverviews; i++ )
510 : {
511 : int nOXSize, nOYSize, iBand;
512 : VWOTInfo *psInfo;
513 : VRTWarpedDataset *poOverviewDS;
514 :
515 : /* -------------------------------------------------------------------- */
516 : /* What size should this overview be. */
517 : /* -------------------------------------------------------------------- */
518 : nOXSize = (GetRasterXSize() + panNewOverviewList[i] - 1)
519 2 : / panNewOverviewList[i];
520 :
521 : nOYSize = (GetRasterYSize() + panNewOverviewList[i] - 1)
522 2 : / panNewOverviewList[i];
523 :
524 : /* -------------------------------------------------------------------- */
525 : /* Create the overview dataset. */
526 : /* -------------------------------------------------------------------- */
527 2 : poOverviewDS = new VRTWarpedDataset( nOXSize, nOYSize );
528 :
529 8 : for( iBand = 0; iBand < GetRasterCount(); iBand++ )
530 : {
531 2 : GDALRasterBand *poOldBand = GetRasterBand(iBand+1);
532 : VRTWarpedRasterBand *poNewBand =
533 : new VRTWarpedRasterBand( poOverviewDS, iBand+1,
534 2 : poOldBand->GetRasterDataType() );
535 :
536 2 : poNewBand->CopyCommonInfoFrom( poOldBand );
537 2 : poOverviewDS->SetBand( iBand+1, poNewBand );
538 : }
539 :
540 2 : nOverviewCount++;
541 : papoOverviews = (VRTWarpedDataset **)
542 2 : CPLRealloc( papoOverviews, sizeof(void*) * nOverviewCount );
543 :
544 2 : papoOverviews[nOverviewCount-1] = poOverviewDS;
545 :
546 : /* -------------------------------------------------------------------- */
547 : /* Prepare update transformation information that will apply */
548 : /* the overview decimation. */
549 : /* -------------------------------------------------------------------- */
550 2 : GDALWarpOptions *psWO = (GDALWarpOptions *) poWarper->GetOptions();
551 2 : psInfo = (VWOTInfo *) CPLCalloc(sizeof(VWOTInfo),1);
552 :
553 2 : strcpy( psInfo->sTI.szSignature, "GTI" );
554 2 : psInfo->sTI.pszClassName = "VRTWarpedOverviewTransform";
555 2 : psInfo->sTI.pfnTransform = VRTWarpedOverviewTransform;
556 2 : psInfo->sTI.pfnCleanup = VRTWarpedOverviewCleanup;
557 2 : psInfo->sTI.pfnSerialize = NULL;
558 :
559 2 : psInfo->pfnBaseTransformer = psWO->pfnTransformer;
560 2 : psInfo->pBaseTransformerArg = psWO->pTransformerArg;
561 :
562 2 : psInfo->dfXOverviewFactor = GetRasterXSize() / (double) nOXSize;
563 2 : psInfo->dfYOverviewFactor = GetRasterYSize() / (double) nOYSize;
564 :
565 : /* -------------------------------------------------------------------- */
566 : /* Initialize the new dataset with adjusted warp options, and */
567 : /* then restore to original condition. */
568 : /* -------------------------------------------------------------------- */
569 2 : psWO->pfnTransformer = VRTWarpedOverviewTransform;
570 2 : psWO->pTransformerArg = psInfo;
571 :
572 2 : poOverviewDS->Initialize( psWO );
573 :
574 2 : psWO->pfnTransformer = psInfo->pfnBaseTransformer;
575 2 : psWO->pTransformerArg = psInfo->pBaseTransformerArg;
576 : }
577 :
578 2 : CPLFree( panNewOverviewList );
579 :
580 : /* -------------------------------------------------------------------- */
581 : /* Progress finished. */
582 : /* -------------------------------------------------------------------- */
583 2 : pfnProgress( 1.0, NULL, pProgressData );
584 :
585 2 : SetNeedsFlush();
586 :
587 2 : return CE_None;
588 : }
589 :
590 : /************************************************************************/
591 : /* GDALInitializeWarpedVRT() */
592 : /************************************************************************/
593 :
594 : /**
595 : * Set warp info on virtual warped dataset.
596 : *
597 : * Initializes all the warping information for a virtual warped dataset.
598 : *
599 : * This method is the same as the C++ method VRTWarpedDataset::Initialize().
600 : *
601 : * @param hDS dataset previously created with the VRT driver, and a
602 : * SUBCLASS of "VRTWarpedDataset".
603 : *
604 : * @param psWO the warp options to apply. Note that ownership of the
605 : * transformation information is taken over by the function though everything
606 : * else remains the property of the caller.
607 : *
608 : * @return CE_None on success or CE_Failure if an error occurs.
609 : */
610 :
611 : CPLErr CPL_STDCALL
612 2 : GDALInitializeWarpedVRT( GDALDatasetH hDS, GDALWarpOptions *psWO )
613 :
614 : {
615 2 : VALIDATE_POINTER1( hDS, "GDALInitializeWarpedVRT", CE_Failure );
616 :
617 2 : return ((VRTWarpedDataset *) hDS)->Initialize( psWO );
618 : }
619 :
620 : /************************************************************************/
621 : /* XMLInit() */
622 : /************************************************************************/
623 :
624 38 : CPLErr VRTWarpedDataset::XMLInit( CPLXMLNode *psTree, const char *pszVRTPath )
625 :
626 : {
627 : CPLErr eErr;
628 :
629 : /* -------------------------------------------------------------------- */
630 : /* Initialize blocksize before calling sub-init so that the */
631 : /* band initializers can get it from the dataset object when */
632 : /* they are created. */
633 : /* -------------------------------------------------------------------- */
634 38 : nBlockXSize = atoi(CPLGetXMLValue(psTree,"BlockXSize","512"));
635 38 : nBlockYSize = atoi(CPLGetXMLValue(psTree,"BlockYSize","128"));
636 :
637 : /* -------------------------------------------------------------------- */
638 : /* Initialize all the general VRT stuff. This will even */
639 : /* create the VRTWarpedRasterBands and initialize them. */
640 : /* -------------------------------------------------------------------- */
641 38 : eErr = VRTDataset::XMLInit( psTree, pszVRTPath );
642 :
643 38 : if( eErr != CE_None )
644 0 : return eErr;
645 :
646 : /* -------------------------------------------------------------------- */
647 : /* Find the GDALWarpOptions XML tree. */
648 : /* -------------------------------------------------------------------- */
649 : CPLXMLNode *psOptionsTree;
650 38 : psOptionsTree = CPLGetXMLNode( psTree, "GDALWarpOptions" );
651 38 : if( psOptionsTree == NULL )
652 : {
653 : CPLError( CE_Failure, CPLE_AppDefined,
654 0 : "Count not find required GDALWarpOptions in XML." );
655 0 : return CE_Failure;
656 : }
657 :
658 : /* -------------------------------------------------------------------- */
659 : /* Adjust the SourceDataset in the warp options to take into */
660 : /* account that it is relative to the VRT if appropriate. */
661 : /* -------------------------------------------------------------------- */
662 : int bRelativeToVRT =
663 : atoi(CPLGetXMLValue(psOptionsTree,
664 38 : "SourceDataset.relativeToVRT", "0" ));
665 :
666 : const char *pszRelativePath = CPLGetXMLValue(psOptionsTree,
667 38 : "SourceDataset", "" );
668 : char *pszAbsolutePath;
669 :
670 38 : if( bRelativeToVRT )
671 : pszAbsolutePath =
672 : CPLStrdup(CPLProjectRelativeFilename( pszVRTPath,
673 36 : pszRelativePath ) );
674 : else
675 2 : pszAbsolutePath = CPLStrdup(pszRelativePath);
676 :
677 38 : CPLSetXMLValue( psOptionsTree, "SourceDataset", pszAbsolutePath );
678 38 : CPLFree( pszAbsolutePath );
679 :
680 : /* -------------------------------------------------------------------- */
681 : /* And instantiate the warp options, and corresponding warp */
682 : /* operation. */
683 : /* -------------------------------------------------------------------- */
684 : GDALWarpOptions *psWO;
685 :
686 38 : psWO = GDALDeserializeWarpOptions( psOptionsTree );
687 38 : if( psWO == NULL )
688 0 : return CE_Failure;
689 :
690 38 : this->eAccess = GA_Update;
691 :
692 38 : if( psWO->hDstDS != NULL )
693 : {
694 0 : GDALClose( psWO->hDstDS );
695 0 : psWO->hDstDS = NULL;
696 : }
697 :
698 38 : psWO->hDstDS = this;
699 :
700 : /* -------------------------------------------------------------------- */
701 : /* Instantiate the warp operation. */
702 : /* -------------------------------------------------------------------- */
703 38 : poWarper = new GDALWarpOperation();
704 :
705 38 : eErr = poWarper->Initialize( psWO );
706 38 : if( eErr != CE_None)
707 : {
708 : /* -------------------------------------------------------------------- */
709 : /* We are responsible for cleaning up the transformer outselves. */
710 : /* -------------------------------------------------------------------- */
711 0 : if( psWO->pTransformerArg != NULL )
712 : {
713 0 : GDALDestroyTransformer( psWO->pTransformerArg );
714 0 : psWO->pTransformerArg = NULL;
715 : }
716 :
717 0 : if( psWO->hSrcDS != NULL )
718 : {
719 0 : GDALClose( psWO->hSrcDS );
720 0 : psWO->hSrcDS = NULL;
721 : }
722 : }
723 :
724 38 : GDALDestroyWarpOptions( psWO );
725 38 : if( eErr != CE_None )
726 : {
727 0 : delete poWarper;
728 0 : poWarper = NULL;
729 : }
730 :
731 : /* -------------------------------------------------------------------- */
732 : /* Generate overviews, if appropriate. */
733 : /* -------------------------------------------------------------------- */
734 : char **papszTokens = CSLTokenizeString(
735 38 : CPLGetXMLValue( psTree, "OverviewList", "" ));
736 : int iOverview;
737 :
738 39 : for( iOverview = 0;
739 : papszTokens != NULL && papszTokens[iOverview] != NULL;
740 : iOverview++ )
741 : {
742 1 : int nOvFactor = atoi(papszTokens[iOverview]);
743 :
744 1 : if (nOvFactor > 0)
745 1 : BuildOverviews( "NEAREST", 1, &nOvFactor, 0, NULL, NULL, NULL );
746 : else
747 : CPLError(CE_Failure, CPLE_AppDefined,
748 0 : "Bad value for overview factor : %s", papszTokens[iOverview]);
749 : }
750 :
751 38 : CSLDestroy( papszTokens );
752 :
753 38 : return eErr;
754 : }
755 :
756 : /************************************************************************/
757 : /* SerializeToXML() */
758 : /************************************************************************/
759 :
760 4 : CPLXMLNode *VRTWarpedDataset::SerializeToXML( const char *pszVRTPath )
761 :
762 : {
763 : CPLXMLNode *psTree;
764 :
765 4 : psTree = VRTDataset::SerializeToXML( pszVRTPath );
766 :
767 4 : if( psTree == NULL )
768 0 : return psTree;
769 :
770 : /* -------------------------------------------------------------------- */
771 : /* Set subclass. */
772 : /* -------------------------------------------------------------------- */
773 : CPLCreateXMLNode(
774 : CPLCreateXMLNode( psTree, CXT_Attribute, "subClass" ),
775 4 : CXT_Text, "VRTWarpedDataset" );
776 :
777 : /* -------------------------------------------------------------------- */
778 : /* Serialize the block size. */
779 : /* -------------------------------------------------------------------- */
780 : CPLCreateXMLElementAndValue( psTree, "BlockXSize",
781 4 : CPLSPrintf( "%d", nBlockXSize ) );
782 : CPLCreateXMLElementAndValue( psTree, "BlockYSize",
783 4 : CPLSPrintf( "%d", nBlockYSize ) );
784 :
785 : /* -------------------------------------------------------------------- */
786 : /* Serialize the overview list. */
787 : /* -------------------------------------------------------------------- */
788 4 : if( nOverviewCount > 0 )
789 : {
790 : char *pszOverviewList;
791 : int iOverview;
792 :
793 1 : pszOverviewList = (char *) CPLMalloc(nOverviewCount*8 + 10);
794 1 : pszOverviewList[0] = '\0';
795 2 : for( iOverview = 0; iOverview < nOverviewCount; iOverview++ )
796 : {
797 : int nOvFactor;
798 :
799 : nOvFactor = (int)
800 : (0.5+GetRasterXSize()
801 1 : / (double) papoOverviews[iOverview]->GetRasterXSize());
802 :
803 : sprintf( pszOverviewList + strlen(pszOverviewList),
804 1 : "%d ", nOvFactor );
805 : }
806 :
807 1 : CPLCreateXMLElementAndValue( psTree, "OverviewList", pszOverviewList );
808 :
809 1 : CPLFree( pszOverviewList );
810 : }
811 :
812 : /* ==================================================================== */
813 : /* Serialize the warp options. */
814 : /* ==================================================================== */
815 : CPLXMLNode *psWOTree;
816 :
817 4 : if( poWarper != NULL )
818 : {
819 : /* -------------------------------------------------------------------- */
820 : /* We reset the destination dataset name so it doesn't get */
821 : /* written out in the serialize warp options. */
822 : /* -------------------------------------------------------------------- */
823 4 : char *pszSavedName = CPLStrdup(GetDescription());
824 4 : SetDescription("");
825 :
826 4 : psWOTree = GDALSerializeWarpOptions( poWarper->GetOptions() );
827 4 : CPLAddXMLChild( psTree, psWOTree );
828 :
829 4 : SetDescription( pszSavedName );
830 4 : CPLFree( pszSavedName );
831 :
832 : /* -------------------------------------------------------------------- */
833 : /* We need to consider making the source dataset relative to */
834 : /* the VRT file if possible. Adjust accordingly. */
835 : /* -------------------------------------------------------------------- */
836 4 : CPLXMLNode *psSDS = CPLGetXMLNode( psWOTree, "SourceDataset" );
837 : int bRelativeToVRT;
838 : char *pszRelativePath;
839 :
840 : pszRelativePath =
841 : CPLStrdup(
842 : CPLExtractRelativePath( pszVRTPath, psSDS->psChild->pszValue,
843 4 : &bRelativeToVRT ) );
844 :
845 4 : CPLFree( psSDS->psChild->pszValue );
846 4 : psSDS->psChild->pszValue = pszRelativePath;
847 :
848 : CPLCreateXMLNode(
849 : CPLCreateXMLNode( psSDS, CXT_Attribute, "relativeToVRT" ),
850 4 : CXT_Text, bRelativeToVRT ? "1" : "0" );
851 : }
852 :
853 4 : return psTree;
854 : }
855 :
856 : /************************************************************************/
857 : /* GetBlockSize() */
858 : /************************************************************************/
859 :
860 73 : void VRTWarpedDataset::GetBlockSize( int *pnBlockXSize, int *pnBlockYSize )
861 :
862 : {
863 73 : assert( NULL != pnBlockXSize );
864 73 : assert( NULL != pnBlockYSize );
865 :
866 73 : *pnBlockXSize = nBlockXSize;
867 73 : *pnBlockYSize = nBlockYSize;
868 73 : }
869 :
870 : /************************************************************************/
871 : /* ProcessBlock() */
872 : /* */
873 : /* Warp a single requested block, and then push each band of */
874 : /* the result into the block cache. */
875 : /************************************************************************/
876 :
877 69 : CPLErr VRTWarpedDataset::ProcessBlock( int iBlockX, int iBlockY )
878 :
879 : {
880 69 : if( poWarper == NULL )
881 0 : return CE_Failure;
882 :
883 69 : const GDALWarpOptions *psWO = poWarper->GetOptions();
884 :
885 : /* -------------------------------------------------------------------- */
886 : /* Allocate block of memory large enough to hold all the bands */
887 : /* for this block. */
888 : /* -------------------------------------------------------------------- */
889 : int iBand;
890 : GByte *pabyDstBuffer;
891 : int nDstBufferSize;
892 69 : int nWordSize = (GDALGetDataTypeSize(psWO->eWorkingDataType) / 8);
893 :
894 : // FIXME? : risk of overflow in multiplication if nBlockXSize or nBlockYSize are very large
895 69 : nDstBufferSize = nBlockXSize * nBlockYSize * psWO->nBandCount * nWordSize;
896 :
897 69 : pabyDstBuffer = (GByte *) VSIMalloc(nDstBufferSize);
898 :
899 69 : if( pabyDstBuffer == NULL )
900 : {
901 : CPLError( CE_Failure, CPLE_OutOfMemory,
902 : "Out of memory allocating %d byte buffer in VRTWarpedDataset::ProcessBlock()",
903 0 : nDstBufferSize );
904 0 : return CE_Failure;
905 : }
906 :
907 69 : memset( pabyDstBuffer, 0, nDstBufferSize );
908 :
909 : /* -------------------------------------------------------------------- */
910 : /* Process INIT_DEST option to initialize the buffer prior to */
911 : /* warping into it. */
912 : /* NOTE:The following code is 99% similar in gdalwarpoperation.cpp and */
913 : /* vrtwarped.cpp. Be careful to keep it in sync ! */
914 : /* -------------------------------------------------------------------- */
915 : const char *pszInitDest = CSLFetchNameValue( psWO->papszWarpOptions,
916 69 : "INIT_DEST" );
917 :
918 69 : if( pszInitDest != NULL && !EQUAL(pszInitDest, "") )
919 : {
920 : char **papszInitValues =
921 19 : CSLTokenizeStringComplex( pszInitDest, ",", FALSE, FALSE );
922 19 : int nInitCount = CSLCount(papszInitValues);
923 :
924 54 : for( iBand = 0; iBand < psWO->nBandCount; iBand++ )
925 : {
926 : double adfInitRealImag[2];
927 : GByte *pBandData;
928 35 : int nBandSize = nBlockXSize * nBlockYSize * nWordSize;
929 35 : const char *pszBandInit = papszInitValues[MIN(iBand,nInitCount-1)];
930 :
931 40 : if( EQUAL(pszBandInit,"NO_DATA")
932 : && psWO->padfDstNoDataReal != NULL )
933 : {
934 5 : adfInitRealImag[0] = psWO->padfDstNoDataReal[iBand];
935 5 : adfInitRealImag[1] = psWO->padfDstNoDataImag[iBand];
936 : }
937 : else
938 : {
939 : CPLStringToComplex( pszBandInit,
940 30 : adfInitRealImag + 0, adfInitRealImag + 1);
941 : }
942 :
943 35 : pBandData = ((GByte *) pabyDstBuffer) + iBand * nBandSize;
944 :
945 35 : if( psWO->eWorkingDataType == GDT_Byte )
946 : memset( pBandData,
947 : MAX(0,MIN(255,(int)adfInitRealImag[0])),
948 18 : nBandSize);
949 32 : else if( adfInitRealImag[0] == 0.0 && adfInitRealImag[1] == 0 )
950 : {
951 15 : memset( pBandData, 0, nBandSize );
952 : }
953 2 : else if( adfInitRealImag[1] == 0.0 )
954 : {
955 : GDALCopyWords( &adfInitRealImag, GDT_Float64, 0,
956 : pBandData,psWO->eWorkingDataType,nWordSize,
957 2 : nBlockXSize * nBlockYSize );
958 : }
959 : else
960 : {
961 : GDALCopyWords( &adfInitRealImag, GDT_CFloat64, 0,
962 : pBandData,psWO->eWorkingDataType,nWordSize,
963 0 : nBlockXSize * nBlockYSize );
964 : }
965 : }
966 :
967 19 : CSLDestroy( papszInitValues );
968 : }
969 :
970 : /* -------------------------------------------------------------------- */
971 : /* Warp into this buffer. */
972 : /* -------------------------------------------------------------------- */
973 : CPLErr eErr;
974 :
975 : eErr =
976 : poWarper->WarpRegionToBuffer(
977 : iBlockX * nBlockXSize, iBlockY * nBlockYSize,
978 : nBlockXSize, nBlockYSize,
979 69 : pabyDstBuffer, psWO->eWorkingDataType );
980 :
981 69 : if( eErr != CE_None )
982 : {
983 0 : VSIFree( pabyDstBuffer );
984 0 : return eErr;
985 : }
986 :
987 : /* -------------------------------------------------------------------- */
988 : /* Copy out into cache blocks for each band. */
989 : /* -------------------------------------------------------------------- */
990 160 : for( iBand = 0; iBand < MIN(nBands, psWO->nBandCount); iBand++ )
991 : {
992 : GDALRasterBand *poBand;
993 : GDALRasterBlock *poBlock;
994 :
995 91 : poBand = GetRasterBand(iBand+1);
996 91 : poBlock = poBand->GetLockedBlockRef( iBlockX, iBlockY, TRUE );
997 :
998 91 : if( poBlock != NULL )
999 : {
1000 91 : if ( poBlock->GetDataRef() != NULL )
1001 : {
1002 : GDALCopyWords( pabyDstBuffer + iBand*nBlockXSize*nBlockYSize*nWordSize,
1003 : psWO->eWorkingDataType, nWordSize,
1004 : poBlock->GetDataRef(),
1005 : poBlock->GetDataType(),
1006 : GDALGetDataTypeSize(poBlock->GetDataType())/8,
1007 91 : nBlockXSize * nBlockYSize );
1008 : }
1009 :
1010 91 : poBlock->DropLock();
1011 : }
1012 : }
1013 :
1014 69 : VSIFree( pabyDstBuffer );
1015 :
1016 69 : return CE_None;
1017 : }
1018 :
1019 : /************************************************************************/
1020 : /* AddBand() */
1021 : /************************************************************************/
1022 :
1023 5 : CPLErr VRTWarpedDataset::AddBand( GDALDataType eType, char **papszOptions )
1024 :
1025 : {
1026 : UNREFERENCED_PARAM( papszOptions );
1027 :
1028 : SetBand( GetRasterCount() + 1,
1029 5 : new VRTWarpedRasterBand( this, GetRasterCount() + 1, eType ) );
1030 :
1031 5 : return CE_None;
1032 : }
1033 :
1034 : /************************************************************************/
1035 : /* ==================================================================== */
1036 : /* VRTWarpedRasterBand */
1037 : /* ==================================================================== */
1038 : /************************************************************************/
1039 :
1040 : /************************************************************************/
1041 : /* VRTWarpedRasterBand() */
1042 : /************************************************************************/
1043 :
1044 : VRTWarpedRasterBand::VRTWarpedRasterBand( GDALDataset *poDS, int nBand,
1045 73 : GDALDataType eType )
1046 :
1047 : {
1048 73 : Initialize( poDS->GetRasterXSize(), poDS->GetRasterYSize() );
1049 :
1050 73 : this->poDS = poDS;
1051 73 : this->nBand = nBand;
1052 73 : this->eAccess = GA_Update;
1053 :
1054 : ((VRTWarpedDataset *) poDS)->GetBlockSize( &nBlockXSize,
1055 73 : &nBlockYSize );
1056 :
1057 73 : if( eType != GDT_Unknown )
1058 7 : this->eDataType = eType;
1059 73 : }
1060 :
1061 : /************************************************************************/
1062 : /* ~VRTWarpedRasterBand() */
1063 : /************************************************************************/
1064 :
1065 73 : VRTWarpedRasterBand::~VRTWarpedRasterBand()
1066 :
1067 : {
1068 73 : FlushCache();
1069 73 : }
1070 :
1071 : /************************************************************************/
1072 : /* IReadBlock() */
1073 : /************************************************************************/
1074 :
1075 : CPLErr VRTWarpedRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
1076 69 : void * pImage )
1077 :
1078 : {
1079 : CPLErr eErr;
1080 69 : VRTWarpedDataset *poWDS = (VRTWarpedDataset *) poDS;
1081 : GDALRasterBlock *poBlock;
1082 :
1083 69 : poBlock = GetLockedBlockRef( nBlockXOff, nBlockYOff, TRUE );
1084 :
1085 69 : eErr = poWDS->ProcessBlock( nBlockXOff, nBlockYOff );
1086 :
1087 69 : if( eErr == CE_None && pImage != poBlock->GetDataRef() )
1088 : {
1089 : int nDataBytes;
1090 : nDataBytes = (GDALGetDataTypeSize(poBlock->GetDataType()) / 8)
1091 0 : * poBlock->GetXSize() * poBlock->GetYSize();
1092 0 : memcpy( pImage, poBlock->GetDataRef(), nDataBytes );
1093 : }
1094 :
1095 69 : poBlock->DropLock();
1096 :
1097 69 : return eErr;
1098 : }
1099 :
1100 : /************************************************************************/
1101 : /* IWriteBlock() */
1102 : /************************************************************************/
1103 :
1104 : CPLErr VRTWarpedRasterBand::IWriteBlock( int nBlockXOff, int nBlockYOff,
1105 8 : void * pImage )
1106 :
1107 : {
1108 : CPLErr eErr;
1109 8 : VRTWarpedDataset *poWDS = (VRTWarpedDataset *) poDS;
1110 :
1111 : /* This is a bit tricky. In the case we are warping a VRTWarpedDataset */
1112 : /* with a destination alpha band, IWriteBlock can be called on that alpha */
1113 : /* band by GDALWarpDstAlphaMasker */
1114 : /* We don't need to do anything since the data will be kept in the block */
1115 : /* cache by VRTWarpedRasterBand::IReadBlock */
1116 8 : if (poWDS->poWarper->GetOptions()->nDstAlphaBand == nBand)
1117 : {
1118 8 : eErr = CE_None;
1119 : }
1120 : else
1121 : {
1122 : /* Otherwise, call the superclass method, that will fail of course */
1123 0 : eErr = VRTRasterBand::IWriteBlock(nBlockXOff, nBlockYOff, pImage);
1124 : }
1125 :
1126 8 : return eErr;
1127 : }
1128 :
1129 : /************************************************************************/
1130 : /* XMLInit() */
1131 : /************************************************************************/
1132 :
1133 : CPLErr VRTWarpedRasterBand::XMLInit( CPLXMLNode * psTree,
1134 66 : const char *pszVRTPath )
1135 :
1136 : {
1137 66 : return VRTRasterBand::XMLInit( psTree, pszVRTPath );
1138 : }
1139 :
1140 : /************************************************************************/
1141 : /* SerializeToXML() */
1142 : /************************************************************************/
1143 :
1144 6 : CPLXMLNode *VRTWarpedRasterBand::SerializeToXML( const char *pszVRTPath )
1145 :
1146 : {
1147 6 : CPLXMLNode *psTree = VRTRasterBand::SerializeToXML( pszVRTPath );
1148 :
1149 : /* -------------------------------------------------------------------- */
1150 : /* Set subclass. */
1151 : /* -------------------------------------------------------------------- */
1152 : CPLCreateXMLNode(
1153 : CPLCreateXMLNode( psTree, CXT_Attribute, "subClass" ),
1154 6 : CXT_Text, "VRTWarpedRasterBand" );
1155 :
1156 6 : return psTree;
1157 : }
1158 :
1159 : /************************************************************************/
1160 : /* GetOverviewCount() */
1161 : /************************************************************************/
1162 :
1163 1 : int VRTWarpedRasterBand::GetOverviewCount()
1164 :
1165 : {
1166 1 : VRTWarpedDataset *poWDS = (VRTWarpedDataset *) poDS;
1167 :
1168 1 : return poWDS->nOverviewCount;
1169 : }
1170 :
1171 : /************************************************************************/
1172 : /* GetOverview() */
1173 : /************************************************************************/
1174 :
1175 1 : GDALRasterBand *VRTWarpedRasterBand::GetOverview( int iOverview )
1176 :
1177 : {
1178 1 : VRTWarpedDataset *poWDS = (VRTWarpedDataset *) poDS;
1179 :
1180 1 : if( iOverview < 0 || iOverview >= poWDS->nOverviewCount )
1181 0 : return NULL;
1182 : else
1183 1 : return poWDS->papoOverviews[iOverview]->GetRasterBand( nBand );
1184 : }
1185 :
|
